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ABSTRACT 


Owing  to  the  lack  of  analytical  and/or  numerical 
tools  with  which  to  properly  investigate  the  resulting 
fiber  orientations  in  a  filled,  low  Reynold's  number  flow 
field,  a  numerical  model  has  been  developed  to  predict 
such  phenomena.  The  motivation  for  such  a  study  is  evi¬ 
dent  from  the  fact  that  the  mechanical  properties  of  a 
molded  two-phase  component  are  inherently  determined  by 
the  fiber  orientations.  Accurately  predicting  the  orien¬ 
tations  of  the  fibers  from  the  rheological  kinematics 
enhances  the  feasibility  of  designing  a  mold  with  the 
aid  of  the  computer. 

Fiber  orientations  in  complex  geometry  flows  of 
dilute  suspensions  are  solved  in  a  two-part  fashion.  The 
initial  step  includes  the  numerical  solution  of  the  steady- 
state  fluid  mechanics  problem  via  the  finite  element  method 
for  either  Newtonian  or  power  law  constitutive  assumptions. 
The  application  of  Jeffery's  theory  [25]  in  conjunction 
with  the  discretized  representation  of  the  velocity  field 
subsequently  determines  the  fiber  orientations.  Numeri¬ 
cally  predicted  fiber  orientations  have  been  verified  with 

iii 


IV 


existing  analytical  solutions  of  a  two-dimensional  channel 
flow  for  Newtonian  and  power  law  fluids.  Application  of 
the  fiber  orientation  model  to  the  expansion  flow  problem 
has  yielded  results  that  correspond  qualitatively  to  those 
found  empirically.  A  complete  set  of  computer  graphics 
routines  has  been  incorporated  with  the  analysis  to  pro¬ 
vide  easy  data  interpretation. 
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1.  PROBLEM  DEFINITION  AND  MOTIVATION 


Industry  is  currently  producing  large  quantities 
of  reinforced  plastics  with  application  in  a  wide  variety 
of  disciplines.  Reinforced  plastics  apply  the  same  prin¬ 
ciple  known  to  the  Egyptians  nearly  three  thousand  years 
ago;  the  fibers  carry  the  load  whereas  the  matrix  acts  as 
a  bonding  material  to  allow  the  fibers  to  maintain  their 
orientation.  By  varying  the  fiber  orientation  one  may  in 
effect  tailor  the  composite  material  in  such  a  manner  to 
suit  its  application.  In  particular,  the  automotive  in¬ 
dustry  views  the  short  fiber  composite  as  replacing  many 
structural  vehicle  components  in  an  effort  to  reduce  weight. 
It  is  a  feasible  solution  to  the  problem  in  light  of  recent 
advances  in  processing  techniques  that  have  shortened  pro¬ 
duction  times.  The  processing  techniques  referred  to  are 
those  appropriate  to  molding  short  fiber  systems;  injection 
molding,  transfer  molding,  and  compression  molding  are  the 
usual  methods  for  producing  components  in  a  reasonable  time, 
so  important  to  high  production  applications. 

In  the  above  processing  techniques,  it  is  assumed 
that  short  fibers  are  combined  with  a  liquid  resin  to  form 
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2 


the  raw  material  (frequently  referred  to  as  the  charge) . 

The  newly  formed  slurry  is  then  forced  to  flow  within  the 
mold  either  through  the  application  of  pressure  or  by 
closing  the  mold.  When  the  material  (thermosets)  has 
filled  the  cavity,  it  is  cured  by  subjecting  it  to  high 
pressure  and  temperature.  Note  that  we  have  chosen  to 
speak  of  both  a  generic  fiber  and  matrix  instead  of  speci- 
fying  the  constituent  phases.  It  is  desired  to  remain  as 
general  as  possible  in  the  problem  formulation  to  accommo¬ 
date  a  broad  class  of  material  systems.  Questions  concern¬ 
ing  specific  material  systems  will  be  undertaken  when  one 
views  the  constitutive  relationships  for  specific  systems. 

The  fiber  orientation  of  composite  materials  is 
known  to  influence  both  modulus  and  strength.  This  fact 
is  especially  evident  in  the  continuous,  layered  systems 
commonly  associated  with  the  aerospace  industry.  However, 
unlike  the  continuous  fiber  systems  where  the  orientation 
is  known,  the  fiber  orientation  of  the  molded  short  fiber 
system  is  yet  to  be  determined.  Its  determination  will 
inherently  depend  upon  the  rheological  behavior  of  the  fi¬ 
brous  suspension.  Thus,  it  seems  that  if  one  is  able  to 
determine  the  relative  position  and  orientation  of  the  fi¬ 
bers,  it  will  be  possible  to  utilize  existing  methodologies 
to  determine:  first,  the  pointwise  material  properties. 
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and  secondly,  the  stress  state  within  the  molded  part. 
Numerical  codes  have  been  developed  to  achieve  the  deter¬ 
mination  of  the  stress  state  once  the  local  material  prop¬ 
erties  are  known.  Without  much  insight,  it  is  easy  to 
see  that  what  is  being  suggested  is  the  possibility  of  a 
computer-aided  mold  design.  The  motivation  is  one  of 
economics,  in  that  the  costly  empirical  development  of 
optimum  part  and  mold  geometry  could  be  eliminated.  No 
longer  would  a  mold  designer  be  obligated  to  fabricate 
the  mold  in  order  to  determine  resulting  fiber  orientation; 
rather  he  would  have  the  ability  to  interact  with  a  mold 
model  and  develop  a  design  accordingly. 

It  was  pointed  out  above  that  both  the  fiber  posi¬ 
tion  and  orientation  are  sought.  Recall  that  questions  of 
material  position  are  usually  connected  with  homogeneity 
while  those  involving  orientation  or  direction  are  asso¬ 
ciated  with  material  anisotropy/isotropy.  It  is  assumed 
that  the  rheology  within  the  mold  will  affect  the  aniso- 
tropy/isotropy  to  a  greater  extent  than  that  of  homoge¬ 
neity.  Most  experimental  data  seem  to  show  that  the  fiber 
distribution  is  relatively  uniform  throughout  a  given 
molded  part,  but  that  the  orientation  is  highly  dependent 
upon  the  flow  conditions.  Thus,  the  major  emphasis  will 
be  placed  upon  determining  the  fiber  orientation  within 
a  given  flow  regime. 


2. 


INTRODUCTION 


The  underlying  motives  for  writing  this  thesis 
are  twofold.  Clearly  the  dominant  reason  is  to  present 
the  findings  from  some  original  research,  as  is  usually 
the  case.  But  secondly,  it  is  important  that  the  manu¬ 
script  be  an  article  of  instruction  as  well.  Care  has 
been  taken  to  insure  that  unfamiliar  ideas  are  developed 
and  useful  examples  are  included. 

It  is  anticipated  that  many  readers  may  not  be 
familiar  with  the  various  methodologies  referred  to  in 
this  text.  For  that  reason,  it  may  be  helpful  to  review 
certain  background  material.  Appendix  A1  provides  an 
introduction  to  fluid  mechanics  which  is  intended  to 
introduce  the  fundamental  concepts  that  are  the  basis 
of  all  work  in  rheology  and  associated  fields.  While 
this  material  may  be  found  in  various  textbooks,  the  con¬ 
cepts  presented  in  Appendix  A1  are  of  a  concise  nature 
and  follow  a  logical  scheme  of  development  pertaining  to 
this  research.  Also  included  in  this  section  are  the 
various  constitutive  equations  that  are  utilized  herein 
as  well.  Several  examples  are  included  to  elucidate  the 
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contrasting  behavior  between  the  Newtonian  and  non- 
Newtonian  fluids. 

Many  of  the  practical  mold  geometries  found  in 
industry  today  are  complex.  This  consideration  has  forced 
the  use  of  a  numerical  scheme  in  order  to  solve  the  fluid 
mechanics  problem.  The  finite  element  method  (FEM)  was 
chosen  primarily  because  of  its  superior  ability  to 
approximate  irregular  boundaries .  The  FEM  as  pertaining 
to  fluid  systems  is  outlined  in  Appendix  A2 ,  making  spe¬ 
cial  mention  of  the  key  considerations,  as  well  as  the 
similarities  that  exist  with  the  solid  mechanics  formu¬ 
lation. 


Intuitive  orientation  models  are  explained  in 
Chapter  3  for  the  shear  and  extensional  flow  fields.  These 
are  especially  useful  in  reconciling  the  orientations  dic¬ 
tated  by  the  analytical  expressions  in  Chapter  4.  In 
Chapter  4  an  expression  for  the  motion  of  a  single  ellip¬ 
soidal  particle  in  a  viscous  flow  is  presented.  This  is 
the  initial  consideration  in  the  overall  process  of  deter¬ 
mining  fiber  orientation  of  a  dilute  suspension  in  an  ar¬ 
bitrary  flow.  Building  upon  this  simplified  model,  one 
extends  the  theory  to  include  discretized  systems  where 
the  flow  field  is  solved  via  the  FEM. 
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Chapter  5  examines  the  numerical  code  COLINS, 
responsible  for  determining  the  fiber  orientation  within 
a  given  mold  from  the  rheological  behavior.  The  detailed 
program  is  examined  to  reveal  its  nature  and  the  conver¬ 
gence  of  the  scheme  is  discussed  briefly.  The  numerical 
code  COLINS  is  then  called  upon  to  solve  several  specific 
examples  in  Chapter  6.  Examples  include  simple  shearing 
flow,  as  well  as  extensional  flows  that  are  found,  say, 
in  a  10:1  expansion.  The  results  of  the  examples  show 
the  mold  geometry  with  individual  fibers  plotted  at  their 
respective  locations  and  orientations.  The  use  of  an 
orientation  parameter  greatly  assists  in  describing  the 
fiber  orientations  in  a  qualitative  manner.  Finally  in 
Chapter  7 ,  several  extensions  to  the  theory  as  motivation 
for  further  research  are  offered. 


3.  QUALITATIVE  ORIENTATION  MODELS 


Included  in  the  discourse  of  Appendices  A1  and  A2 , 
material  has  been  provided  to  construct  the  necessary  back¬ 
ground  needed  to  appreciate  the  objective  of  this  work.  At 
this  point  it  is  our  intent  to  promulgate  some  of  the  intu¬ 
itive  ideas  concerning  the  orientation  of  fiber-like  parti¬ 
cles  in  a  fluid  flow.  The  models  which  are  discussed  are 
based  on  the  general  trends  exhibited  by  suspended  parti¬ 
cles  that  have  been  observed  experimentally.  This  quali¬ 
tative  appeal  toward  understanding  the  mechanisms  will  be 
useful  later,  when  one  endeavors  to  analyze  fiber  orien¬ 
tation  behavior  for  geometrically  complicated  domains. 


The  last  few  years  have  seen  numerous  experimental 
studies  that  have  illustrated  the  orientation  trends  of 
fiber-filled  systems  undergoing  characteristic  flow  patterns. 
The  first  such  orientation  model  exists  when  the  suspension 
is  subject  to  a  shearing- type  motion.  Regions  of  shearing 
flow  are  characterized  by  a  velocity  gradient  tensor  of 
the  following  form: 
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where  F  is  defined  as  the  shear  rate.  The  coordinate 
system  that  has  been  adopted  for  this  discussion  finds 
the  x1  axis  parallel  to  the  flow  direction  while  the  x2 
axis  is  aligned  with  the  direction  of  varying  velocity. 

For  completeness,  the  x^  axis  is  so  arranged  to  form  a 
right-handed  system  of  the  three  coordinate  directions. 
Under  steady,  two-dimensional  flow,  the  deformation  rate 
may  depend  upon  the  spacial  coordinate,  (x2) ,  but  not 
upon  time . 

Common  occurrences  of  shearing  flow  are  found  in 
the  vicinity  of  physical  boundaries  such  as  mold  walls 
or  inserts.  Clearly  this  is  so  in  light  of  the  "no  slip" 
boundary  condition  at  the  wall,  which  requires  the  local 
velocity  of  the  fluid  to  be  zero  at  the  wall  surface  and 
yet  attain  a  finite  value  in  the  bulk  of  the  flow.  Figure 
3.1  depicts  a  single  fiber,  initially  transverse  to  the 
flow  direction,  in  a  shearing  flow.  if  the  local  velocity 


fiber 


Figure  3.1  Orientation  from  shearing  flow 
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field  surrounding  the  fiber  is  separated  into  two  parts, 
it  is  seen  from  the  figure  that  the  distinct  contribu¬ 
tions  may  be  represented  as  translational  and  rotational 
velocity  components,  respectively.  The  translational 
velocity  acts  only  to  transport  the  fiber  with  the  mean, 
local  fluid  velocity.  Thus,  since  the  particle  center 
moves  with  the  same  velocity  as  the  fluid  there  can  be 
no  resultant  force  acting  on  it.  The  only  possible  action 
on  the  particle  results  from  a  couple,  which  consists  of 
two  contributions;  the  first  tends  to  rotate  the  fiber 
with  the  same  angular  velocity  of  the  surrounding  fluid 
while  the  second  acts  to  align  the  particle  axes  with 
those  of  the  principal  axes  of  the  deformation  rate  ten¬ 
sor.  In  view  of  these  considerations  it  is  shown  that  the 
equilibrium  position  for  a  large  aspect  ratio  ellipsoid 
is  parallel  to  the  fluid  streamlines.  Figure  3.2  shows 
experimentally  the  resultant  fiber  orientation  in  regions 
of  steady  shearing  flow. 

Equally  important  to  realize  is  the  mechanism 
responsible  for  orienting  the  fibers  where  the  flow  pos¬ 
sesses  significant  dilational  components  of  the  deforma¬ 
tion  rate  tensor.  For  convenience  in  describing  flows 
of  this  type,  the  radial  coordinate  system  (r,0)  as  shown 
in  Figure  3.3  is  introduced.  In  terms  of  this  coordinate 
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direct  ion  ^ 
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end  gate 


Figure  3.2  Experimentally  observed  fiber  orientation 
in  shear  flow  for  the  P.E.T.  (polyethylene 
terephthalate)  material  system  (100x) 
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(a)  (b) 

Figure  3.3  Orientation  from  elongation 


system,  then,  the  deformation  rate  tensor  appears  as 

11  0 
o  d22  0 

0  0  0 

and  the  corresponding  stresses  acting  on  the  fluid  ele¬ 
ment  are  also  shown  in  Figure  3.3.  The  radially  directed 
stress  decelerates  the  fluid  in  that  coordinate  direction, 
and  from  the  equations  of  motion,  the  tangential  compo¬ 
nent  of  stress  tends  to  elongate  a  fluid  element  in  the 
0-direction.  If  one  imagines  that  a  fiberlike  particle 
is  within  a  fluid  element  undergoing  such  a  de¬ 
formation,  then  it  is  not  hard  to  visualize  the  appro¬ 
priate  motion  of  the  particle.  Figure  3.3  depicts  the 
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SQui librium  position  of  a  fiber  subjected  to  an  elongated 
flow  situation.  Experimentally  this  effect  has  been  ob¬ 
served  in  the  molding  of  numerous  short  fiber  composites, 
^isure  3.4  illustrates  orientation  patterns  as  the  slurry 
enters  a  region  of  expanding  flow.  Clearly  the  flow  in 
this  region  approximates  the  flow  situation  of  Figure  3.3 
and  as  a  result  tangential  orientation  is  anticipated. 

In  contrast  with  the  equilibrium  position  for 
expansion  flows  described  above ,  converging  flow  will 
orient  a  particle  parallel  to  fluid  streamlines.  To  see 
this,  one  notes  the  stresses  on  a  fluid  element  act  in 
opposite  directions  to  those  of  Figure  3.3.  The  intui¬ 
tive  arguments  for  radially  oriented  fibers  in  converging 
flow  are  outlined  in  Figure  3.5. 

So  far,  we  have  alluded  to  the  equilibrium  posi¬ 
tion  of  fibers  undergoing  certain  flow  criteria.  Nothing 
has  been  mentioned  concerning  the  transient  response  of 
the  suspended  particles —  i.e.  how  do  they  reach  their 
equilibrium  state.  The  next  chapter  will  be  concerned 
with  investigating  the  motion  of  the  particles  from  a 
quantitative  point  of  view. 


13 


Figure  3.4  Experimentally  observed  fiber 

orientation  in  region  of  diverg¬ 
ing  flow  (100x) 


Figure  3 . 5  Orientation  model  for  convering 
flow 


4.  QUANTITATIVE  ORIENTATION  MODELS 


4 . 1  Jeffery's  Orientation  Equation 

To  this  point  one  will  have  noted  certain  quali¬ 
tative  trends  in  particle  orientation  that  have  been  ob¬ 
served  experimentally.  Explanations  of  these  trends  were 
offered  in  the  previous  section.  Within  this  section  the 
basic  equations  and  assumptions  needed  to  describe  the 
particle  motions  quantitatively  are  introduced. 

The  earliest  work  addressing  the  problem  of  sus¬ 
pended  particle  motions  was  by  G.  B.  Jeffery  [25].  Here, 
Jeffery  endeavored  to  solve  the  linearized  Newtonian 
equations  of  motion  for  a  single  ellipsoidal  particle 
immersed  in  a  viscous  flow.  The  boundary  conditions  he 
proposed  were  the  following:  i)  the  fluid  shall  adopt 
the  same  motion  as  that  of  the  particle  at  the  solid- 
fluid  interface  and  the  fluid  shall  attain  the  free  stream 
velocity  at  sufficient  distances  away  from  the  particle; 
ii)  the  free  stream  velocity  was  taken  to  be  a  general 
shearing-type  flow  with  a  fully  populated  deformation 
rate  tensor.  In  his  early  development  of  the  problem. 
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Jeffery  assumed  that  the  components  of  the  deformation 
rate  tensor  were  constant  in  space  on  a  scale  that  is 
large  when  compared  with  the  particle  dimensions.  Later 
this  restriction  is  modified  to  include  flows  such  that 
the  components  of  the  deformation  rate  tensor  are  in¬ 
deed  functions  of  the  spacial  coordinates.  The  orig¬ 
inal  solution  is  carried  out  formally  in  an  ellipsoidal 
coordinate  system.  Here  the  axes  of  the  ellipsoid  serve 
to  define  a  local,  orthogonal  coordinate  system  x^°, 
x2° '  x3° '  trans^-a't-'-nJ  and  rotating  with  the  particle, 
such  that  x ^/b ^  +  x^/b^  +  x°^/b^  =  1  .  Rather  than 

discussing  the  solution  methods  here,  the  reader  is 
referred  to  the  original  paper  for  details  [25] .  The 
results  of  the  analysis  assert  that  the  resultant  force 
on  the  particle  vanishes  as  it  translates  with  the  local, 
mean  velocity  of  the  fluid.  Secondly,  the  resultant 
couple  acting  on  the  particle  consists  of  two  parts: 
one  part  is  due  to  the  local  rotation  (vorticity)  of  the 
fluid  and  the  other  results  from  the  fluid  deformation. 

It  is  natural  to  define  the  components  of  the 
resultant  torque  in  terms  of  a  fixed  or  inertial  coor¬ 
dinate  system,  since  the  fluid  kinematics  will  most 
naturally  be  defined  in  the  latter  reference  frame. 

For  this  purpose,  then,  the  Euler  angles  and  0^  , 
defined  by  Figure  4.1,  are  introduced. 
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v3=rx2 


Figure  4.1  Coordinate  system  used  in 
describing  fiber  rotation 


From  the  expressions  denoting  the  resultant  couple, 
it  may  then  be  shown  [18]  that  the  fiber  motion  may  be 
defined  in  terms  of  the  Euler  angles  such  that 


Z1  ~  Z2cos<})1cotei  -  Z^sincfi^cote-^ 

+  B  [-d^sin^cote^  +  d^cosC^ 

+  d^costf^cotS^  -  y2(d22-d  )sin2<|>. 


(4.1) 


-Z2sin<}>^  +  Z^coscj)^  +  B  [d12cosd^cos20^ 

+  1/2d23sin2(f)xsin20l  + 

+  14(d22-d33)cos2^1sin201 
+  34  (d22+d33)  sin291] 


(4.2) 
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where  the  dot  ( • )  denotes  differentiation  with  respect 
to  time.  Also 


=  components  of  the  vorticity  vector 

[zi  =  1/2(uj>k‘uk,j)eijkI 

=  components  of  the  deformation  rate  tensor 

[dij  “  1/2  (ui,  j+uj.i)  1 

=  (r  2  -  l)/(r  2  +  1) 

P  P 

=  b^/b2  defines  the  particle  aspect  ratio. 


If  one  confines  the  analysis  to  two-dimensional 
fluid  flows,  such  that  variations  in  the  third  coordinate 
direction  may  be  ignored,  then  the  vorticity  vector  will 
be  perpendicular  to  the  plane  of  the  flow.  In  light  of 
this  limitation  eq.(4.1)  may  be  simplified  to 

+  B  [d23cos2(|>1-  1/2(d22-d33)sin2(^1]  (4.3) 

Eq.(4.2)  is  no  longer  needed  to  describe  the  in-plane 
fiber  orientation  thus  it  shall  be  disregarded  in  further 
discussion . 


To  reiterate,  eq.(4.3)  describes  the  in-plane 
orientation  response  of  a  single-ellipsoidal  particle 
in  a  general  Newtonian  shearing  flow.  Within  the  pre¬ 
vious  statement  lie  the  bulk  of  the  assumptions  and 
limitations  of  the  quantitative  descriptive  ability  of 
eq.C4.3).  Even  in  its  simplified  form  eq.(4.3)  is 
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difficult  to  solve  except  for  situations  where  the 
vorticity  and  deformation  rates  are  easily  expressible 
in  terms  of  the  flow  field.  The  following  sections 
provide  examples  of  two  such  flow  fields  that  allow 
analytical  expressions  to  be  found  for  the  components 
of  d  and  .  In  these  particular  cases,  eq.(4.3)  may 
be  integrated  directly  to  yield  analytical  expressions 
for  the  fiber  orientation  as  a  function  of  time. 

4 . 2  Solution  of  the  Orientation  Equation  by 

Analytical  Methods 

4.2.1  2-D  Poiseuille  Flow 

Consider  a  pressure-driven,  wall-bounded  shear 
flow  to  exist  such  that  the  velocity  variation  occurs 
in  the  X2-x3  Plane*  Figure  4.2  illustrates  the  nomen¬ 
clature  that  will  be  utilized  for  the  subsequent  dis¬ 
cussion.  x. 


Figure  4.2  Two-dimensional  wall-bounded 
shearing  flow 
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For  a  unit  flow,  the  normalized  velocity  profile  has 
been  shown  in  eq.(A1.65)  to  be 


a 


(4.4) 


and  the  components  of  the  vorticity  vector  and  deforma¬ 
tion  rate  tensor  are  respectively  calculated  to  be 


Z.  = 

1 


0 


0 


(4.5) 


0  0  0 
o  o  y2r 
o  y2r  o 


(4.6) 


where  T  =  =  -3x„/a2  . 

ax2  z 

Inserting  eqs.(4.5)  and  (4.6)  into  eq.(4.3),  the 
expression  for  the  angular  orientation  of  the  particle 
becomes 

0,  =  — y —  (r  2  cos2<f>  +  sin2^.)  .  (4.7) 

x  r  +1  p  x  x 

P 

Recall  that  r^  was  defined  as  the  aspect  ratio  of  the 
particle  (b^/b2) .  One  nominally  chooses  rp = 50  for  the 
discussion  of  fiber-like  particles.  Integration  of 
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eq.  (4.7)  with  respect  to  time  admits  the  solution 


tand).  =  r  tan 
1  P 


rt 


r  +  — 

P 


(4.8) 


If  one  designates  that  the  initial  fiber  orientation 
is  perpendicular  to  the  flow  (i.e.  ^  =  0  ,  t=0)  then 

the  corresponding  particle  response  is  shown  in  Figure 
4.3.  The  abscissa  is  a  product  of  the  shear  rate  and 
the  time;  for  a  constant  shear  rate  Figure  4.3  shows 
the  variation  of  <j>^  as  a  function  of  time.  However, 
for  a  non-uniform  shear  rate  some  explanation  is  needed 
to  interpret  the  particle  response  in  flows  of  this  type. 
In  general,  particles  will  assume  orientations  that  de¬ 
pend  upon  their  relative  position  in  the  flow,  as  well 
as  the  duration  of  time. 


First  and  foremost,  it  was  previously  pointed  out 
that  Jeffery's  original  assumptions  included  a  deforma¬ 
tion  rate  tensor  whose  components  were  spacially  inde¬ 
pendent.  Now  this  restriction  is  removed  in  light  of  the 
work  of  Goldsmith  and  Mason  (18]  who  showed  quantitative 
agreement  with  the  rotational  motions  of  particles  between 
Couette  flow  and  Poiseuille  flow.  Provided  that  the  par¬ 
ticles  are  small  when  compared  to  the  channel  half-width, 
the  rotations  proved  identical  between  the  two  flows 


despite  the  linear  variation  of  shear  rate  for  Poiseuille 
flow.  Hence,  Figure  4.3  may  also  be  interpreted  as  pro¬ 
viding  high  degrees  of  orientation  for  particles  subjected 
to  high  shear  rates  for  short  durations  of  time,  as  it  is 
the  product  ft  that  determines  the  orientation. 


i  5  10 

/r/t 

Figure  4.3  Angular  orientation  of  a  particle 
subjected  to  shearing  flow 
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If  one  considers  an  alternative  method  of  inter¬ 
preting  Figure  4.3  we  define  the  quantity 

D  =  | T | t  (4.9) 


that  may  be  regarded  as  a  measure  of  the  degree  of 
equilibrium  particle  orientation  for  shearing  flows. 

We  shall  endeavor  to  illustrate  the  behavior  of  D  through¬ 
out  a  flow  regime  to  isolate  the  regions  of  various  de¬ 
grees  of  orientation.  Toward  this  end,  it  is  advanta¬ 
geous  to  recast  the  right-hand  side  of  eq.(4.9)  in  terms 
of  the  coordinate  positions  within  the  domain.  Suppose 
that  the  local  velocity  is  given  by 
dx^  x 

u  =  =  —  (for  steady  flow)  (4.10a) 

X3 

t  =  —  (4.10b) 


where  one  disallows  velocity  variations  with  x^  .  Utiliz¬ 
ing  eq. (4.10b)  and  substituting  for  the  local  velocity 
from  eq.(4.4),  the  expression  for  D  (eq.(4.9))  becomes 


3/2  [1  ~  (x22/a2)  ] 

or  alternatively, 

D (a2-x  2) 

x_  =  - =■ - 

3  2x2 


2x2X3 


(4.11) 


(4.12) 
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For  convenience,  the  dimensionless  channel  coor¬ 
dinates  are  introduced  in  the  form  of  eqs. (4.13a)  and 
(4.13b) 

*3  "IT  (4.13a) 

*2  'IT  (4.13b) 

When  these  modifications  are  added  to  eq.(4.12)  the 
resulting  expression  becomes 
D(l-x  2) 

x,  =  - —  (4.14) 

3  2x2 

The  consequence  of  eqs. (4.9)  to  (4.14)  has  been  to 
remove  the  time  dependence  from  the  orientation  equa¬ 
tion  (eq.(4.8))  in  exchange  for  a  spacial  dependence. 

This  was  accomplished  explicitly  in  eq. (4.10b).  If 
one  plots  the  behavior  of  x^  as  a  function  of  x2  with 
D  as  a  parameter,  the  results  are  given  in  Figure  4.4a. 

To  interpret  the  figure,  the  reader  is  reminded  that 
the  value  D  is  a  measure  of  the  degree  of  equilibrium 
orientation  (see  Figure  4.3).  The  flow  regime  consists 
of  half  a  channel  with  an  aspect  ratio  (x2/x2)  of  10  . 
Regions  bounded  by  contour  lines  indicate  the  respective 
areas  of  orientation  that  result  from  the  flow  conditions. 
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flow  direction 


Figure  4.4a  Variation  of  the  degree  of  equilibrium 
orientation  for  2-D  Poiseuille  flow 


Figure  4.4b  Three-dimensional  plot  of  D=f(x„,x,) 
for  2-D  Poiseuille  flow  (n  =  l)  J 
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Another  way  to  summarize  the  above  discussion  is 
to  recognize  that  eq.(4.14)  may  be  solved  explicitly 
for  D  . 


D 


2x2X3 


(1-x ,) 


(4.15) 


From  eq.(4.15)  one  may  note  that  D  is  a  function  of  both 
x^  ana  x^  ,  and  may  by  represented  as  a  surface  in  three- 
dimensional  space.  This  surface  is  plotted  in  perspective 
in  Figure  4.4b,  with  the  variation  of  D  being  plotted  out 
of  the  x^ -x^  plane.  From  eq.(4.15),  the  singular  behav¬ 
ior  of  the  function  at  x2  =  1  may  easily  be  seen.  Factor¬ 
ing  the  denominator,  the  strength  of  the  singularity  is 
recognized  to  be  of  the  first  order  for  the  case  of  New¬ 
tonian  flow.  The  singularity  at  the  wall  (x^  =  1)  may  be 
attributed  to  the  fact  that  the  velocity  vanishes  there 
resulting  in  a  division  by  zero  in  eq. (4.11) .  The  extreme 
values  of  the  parameter  D  near  the  wall  reiterate  the 
notion  of  highly  aligned  fibers  parallel  to  the  flow  ad¬ 
jacent  to  the  boundary  of  the  domain. 


4.2.2  Power  Law  Flow 

If  one  considers  the  same  domain  as  shown  in 
Figure  4.2,  but  changes  the  constitutive  behavior  of  the 
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fluid  in  the  channel,  the  effect  of  a  shear  thinning 
(pseudoplastic)  flow  on  particle  behavior  can  be  inves¬ 
tigated.  The  flow  remains  steady  in  the  x^  direction 
and  possesses  a  nonuniform  velocity  gradient  in  the  x 2 
direction  (see  Figure  Al.ll  for  n<l) .  This  problem 
is  governed  by  the  same  assumptions  of  the  previous 
section,  where  nonuniform  velocity  gradients  were  shown 
to  correctly  predict  the  particle  orientation. 


The  normalized  velocity  profile  for  a  unit  flow 
of  a  power  law  fluid  can  be  shown  to  be  (see  Appendix 


Al.  3) 


u 


1  +  2n 
1  +  n 


[1  -  (x2/a) 


n+1 


(4.16) 


where  n  =  power  law  index.  Straightforward  differention 
yields  the  particular  component  of  the  deformation  rate 
tensor,  namely 


n 

(1 +  2n)  X2 

n  n+1 

n 

a 


(4.17) 


Following  the  identical  sequence  of  events  of  the  previous 
section,  one  finds  the  degree  of  orientation  to  be 


D 


r  1 1 


(1  +  n) 
n 


1 

,  n  1 

X3X2  n+1  n+1 


(4.18) 
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Introducing  the  normalized  coordinates  of  eqs.(4.13), 

the  desired  expression  analogous  to  eq.(4.14)  becomes 

fn+l) 


x. 


=  D 


n 


1  +  n 


(1  -  x. 


n 


(4.19) 


Since  it  is  of  interest  to  discern  the  effect  of  the 
pseudoplasticity  of  the  fluid  on  particle  orientations, 
the  value  of  D  is  taken  as  constant  and  n  is  the  para¬ 
meter.  In  Figure  4.5a,  a  plot  of  eq. (4.19)  is  shown 
where  the  effect  of  the  power  law  constitutive  equation 
may  be  seen  upon  the  fiber  orientation.  This  result  can 
be  explained  if  one  recalls  the  associated  velocity  pro¬ 
file  for  shear  thinning  flows  in  a  channel  (see  Figure 
Al.ll).  For  these  types  of  flow  conditions  there  exists 
a  region  where  the  velocity  profile  is  flat,  originating 
at  the  tube  centerline  and  extending  radially  outward. 
Corresponding  to  this  flat  velocity  profile  region  is 
a  random  core  of  fibers.  Oriented  regions  of  fibers  are 
found  only  in  a  boundary  layer  near  the  wall  whose  thick¬ 
ness  varies  according  to  the  index  n  .  These  distributions 
of  fibers  are  characteristic  of  shear  thinning  flows  and 
have  been  observed  experimentally. 


For  purposes  of  comparison.  Figure  4.5b  illus¬ 
trates  the  continuous  variation  of  D  throughout  the  x 2  , 
x^  domain  of  Figure  4.5a.  The  discrete  value  for  the 
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Figure  4.5b  Three-dimensional  plot  of  D 
for  power  law  fluid  in  wall 
shear  flow  (n=0.1) 
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power  law  index  was  chosen  to  be  0.1,  which  approxi¬ 
mates  a  highly  shear  thinning  flow.  The  important 
comparison  to  note  arises  from  the  variation  of  D 
close  to  the  wall.  Here,  it  may  be  shown  that  for  a 
general  power  law  fluid,  the  singularity  occurring  at 
x2  =  1  is  of  the  order 


0(1/(1  -  x2 


n+1 


2n 


(4.20) 


This  gives  rise  to  a  distinct  boundary  layer  and  identi¬ 
fiable  core  region. 


So  far,  it  has  been  assumed  that  the  initial 
fiber  orientations  have  been  perpendicular  to  the  fluid 
streamlines.  Suppose  that  the  initial  fiber  distribu¬ 
tion  is  assumed  to  be  random  in  nature.  What  will  the 
resulting  fiber  distribution  be  within  the  half  channel? 

A  plot  of  this  situation  is  shown  in  Figures  4.6a -d  , 
for  decreasing  values  of  the  power  law  index.  Clearly, 
these  illustrations  of  fiber  orientation  display  the 
aligned  regions  of  fibers  and  enhance  the  notion  of  a 
distinct  boundary  layer  for  small  values  of  n  .  This  close 
agreement  among  Figures  4.5a  and  4.6a -d  indicates  that 
the  initial  conditions  have  a  negligible  effect  on  the 
downstream  part  of  the  solution.  It  is  the  flow  kinema¬ 
tics  that  exercise  a  major  role  in  the  resulting  orienta¬ 
tion  of  the  particles  for  these  dilute  systems. 


(Figure  4.6  continued) 
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4 . 3  Application  to  the  Discretized  System 

To  this  point  simplified  flows  have  been  dis¬ 
cussed  such  that  the  velocity  fields  (or  derivatives 
thereof)  have  been  expressed  by  analytical  functions. 

This  was  necessary  to  facilitate  the  closed-form  inte¬ 
gration  of  eq.(4.3) .  We  wish  to  extend  the  ideas  pre¬ 
sented  thus  far  to  more  complicated  domains,  where  the 
velocity  fields  may  not  be  solved  exactly.  The  exten¬ 
sion  is  obvious  as  one  seeks  to  model  the  orientation 
patterns  that  exist  in  complex  molds . 

The  flow  chart  of  Figure  4.7  illustrates  how  the 
method  of  solution  is  altered  when  the  flow  becomes  com¬ 
plex  due  to  geometry.  The  initial  difference  occurs  upon 
implementation  of  the  FEM  to  solve  the  fluid  mechanics 
problem.  As  a  result  of  this  approximation  to  the  velo¬ 
city  field,  one  maintains  only  a  discretized  representa¬ 
tion  of  components  of  velocity.  Derived  quantities  such 
as  stress  components,  stream  function,  and  vorticity  will 
also  be  of  a  discretized  nature;  thus,  all  quantities  of 
interest  are  only  known  at  the  nodal  points  within  the 
domain. 


Since  the  extension  of  this  analysis  is  to  arbitrary 
flows,  one  may  not  assume  a  priori  the  components  of  the 
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Figure  4.7  Alternative  methods  of  determining  fiber 
orientations  in  dilute  systems 
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deformation  rate  tensor  to  be  constant  in  space  on  a 
scale  that  is  large  in  comparison  with  the  particle. 
However,  by  resorting  to  the  discretized  domain  the 
original  problem  of  determining  fiber  orientation  has 
been  recast  into  a  series  of  independent,  smaller  prob¬ 
lems.  For  each  subproblem  the  components  of  distortion 
and  rotation  are  taken  to  be  constant  in  space  and  the 
original  assumptions  of  Jeffery  are  recovered.  The 
numerical  integration  of  the  fiber  orientation  equation 
(eq.(4.3))  along  a  streamline  admits  a  piecewise  contin¬ 
uous  approximation  to  the  true  solution.  The  components 
of  the  deformation  rate  tensor  and  vorticity  vector  are 
updated  at  each  stage  of  the  integration  but  remain  con¬ 
stant  during  the  integration.  The  convergence  of  such 
a  scheme  is  investigated  in  section  5.2. 

In  the  following  chapter,  the  methodology  of 
solving  for  particle  orientations  when  only  a  discre¬ 
tized  velocity  field  is  known  is  presented  in  some  de¬ 
tail.  As  in  other  sections,  the  subsequent  chapter  draws 
on  important  discussions  in  the  appendices  and  will  cul¬ 
minate  with  the  ability  to  determine  the  fiber  orienta¬ 
tions  within  molds  of  complicated  designs. 


5.  THE  NUMERICAL  ROUTINE  "COLINS 


5 . 1  Methodology 

The  emphasis  in  this  chapter  will  be  placed  upon 
examining  the  numerical  program  COLINS  in  detail.  Incor¬ 
porated  within  the  program  is  the  methodology  allowing 
one  to  predict  and  view  fiber  orientations  within  a  gen¬ 
eral  flow  system.  It  is  assumed  that  the  fluid  mechanics 
problem  has  been  completely  solved  and  a  discrete  repre¬ 
sentation  of  the  primitive  and  derived  quantities  is  in 
hand.  Thus,  the  input  to  the  program  COLINS  is  a  com¬ 
plete  knowledge  of  the  flow  problem.  For  the  present 
analysis,  it  is  assumed  that  the  fiber  orientations  will 
not  alter  the  velocity  field.  This  criterion  may  be  sat¬ 
isfied  if  one  restricts  his  considerations  to  dilute  sus¬ 
pensions.  Usually  the  standard  for  dilute  systems  is 
taken  to  be  no  greater  than  0.05  particulate  volume  frac¬ 
tion.  To  correctly  model  the  resulting  fiber  orientations 
in  a  more  concentrated  material  system,  an  iterative-type 
solution  is  required.  Here,  the  effects  of  fiber  orien¬ 
tation  on  viscosity  are  predicted  and  the  velocity  field 
is  recalculated  in  light  of  this  updated  viscosity  function. 
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Returning  to  the  problem  at  hand,  one  seeks  to 
provide  a  means  of  solving  eq.(4.3)  within  the  discre¬ 
tized  domain.  The  first  consideration  assumes  that  the 
center  of  gravity  of  the  fiber  will  translate  parallel 
to  fluid  streamlines.  Streamlines  are  easily  calculated 
from  the  discretized  representation  of  the  stream  func¬ 
tion  as  follows.  The  entire  grid  is  scanned  element  by 
element  to  identify  the  nodal  points  corresponding  to  a 
given  value  of  the  stream  function  (tp)  .  More  than  likely 
the  desired  value  of  ip  will  not  intersect  an  element 
boundary  precisely  at  a  nodal  point.  Linear  interpolation 
is  used  to  pinpoint  the  precise  location,  and  hence  the 
coordinate  values,  for  the  desired  value  of  ip .  Upon  scan¬ 
ning  the  entire  array  of  elements,  a  collection  of  x,y 
coordinate  values  is  obtained  that  defines  a  particular 
streamline.  The  next  task  requires  that  these  points  be 
organized  appropriately  so  as  to  define  a  path  line  of  a 
particle  in  the  direction  of  the  flow. 

Scanning  element  by  element  will  duplicate  all  of 
the  stream  function  points  excluding  the  first  one  and 
the  last.  The  reason  for  this  can  be  seen  in  Figure  5.1, 
where  it  is  noticed  that  if  a  streamline  intersects  an 
element,  it  must  intersect  an  element  boundary  twice. 

Given  the  initial  intersecting  point  of  a  streamline,  Pq  , 
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Figure  5.1  Constructing  a  path  line 

and  its  neighboring  point,  ,  the  subroutine  ARRANG 
(of  the  program  COLINS)  searches  until  it  locates  P^'  and 
hence  its  neighbor  P ^  •  The  procedure  continues  until 
an  entire  path  line  is  constructed  along  which  a  parti¬ 
cle  will  translate. 

Now  that  a  typical  particle  path  is  known,  eq. 

(4.3)  may  be  numerically  integrated  along  the  path  line. 

Since  the  independent  variable  is  time  in  eq.(4.3),  we 

need  to  interpret  the  distance  along  the  streamline  in 

terms  of  time  steps.  Refer  to  Figure  5.1  and  consider 

the  distance  between  ip  and  ip.  .  If  one  divides  the  dis- 

o  1 

tance  by  the  local  fluid  velocity  at  f  ,  a  natural  time 
step  is  defined.  Clearly  the  accuracy  of  this  procedure 
is  determined  by  the  fineness  of  the  mesh  and  is  inves¬ 
tigated  in  the  next  section. 
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The  components  of  the  vorticity  and  distortion 
tensors  have  been  assumed  to  vary  with  position  in  an 
arbitrary  flow.  More  particularly,  they  will  vary  along 
the  streamline  during  integration.  This  situation  may 
be  easily  handled  if  the  values  of  and  d  are  updated 
at  each  time  step  in  the  integration.  Again,  the  pre¬ 
cise  values  for  and  d  at  the  points  ifm  (Figure  5.1) 
are  determined  via  linear  interpolation  between  the  nodal 
values.  Hence,  the  fiber  orientation  may  be  numerically 
integrated  via  a  canned,  first-order,  nonlinear  equation 
solver.  Fiber  orientations  may  be  solved  for  along  many 
streamlines  which  results  in  the  fiber  orientations  being 
known  for  the  entire  flow  field. 

5 . 2  Convergence  of  the  Numerical  Scheme 

In  this  section,  the  question  of  convergence  of  the 
numerically  integrated  orientation  equation  (eq.(4.3))  is 
investigated.  Attention  is  focused  entirely  upon  the  above 
questions,  thus  assuming  a  convergent  fluid  flow  field  al¬ 
ready  exists.1  Convergence  in  the  context  of  linear  fluid 
mechanics  problems  is  defined  as  approaching  the  exact 

1This  assumption  will  be  satisfied  by  considering  two- 
dimensional  Poiseuille  flow,  a  problem  whose  finite- 
element  solution  yields  the  exact  analytical  represen¬ 
tation  of  the  velocity  and  pressure  fields. 
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solution  by  successive  mesh  refinements.  This  idea  of 
convergence  will  be  adopted  here,  to  treat  the  problem 
at  hand. 


Recall  that  the  fiber  orientation  is  solved  for 
by  numerically  integrating  eq.(4.3)  along  selected  path 
lines.  The  coefficients  of  this  first-order,  nonlinear 
equation  are  updated  at  each  point  where  the  path  line 
intersects  an  element  boundary.  It  stands  to  reason  that 
the  exact  solution  will  be  more  closely  approximated  if 
there  is  an  increased  number  of  such  points.  Hence,  it 
follows  that  a  continuous  refinement  of  the  mesh  should 
lead  to  convergence  of  the  fiber  orientation  solution. 

Figures  5.2a-  c  show  the  three  respective  grids 
for  2-D ,  symmetric  Poiseuille  flow.  The  discretized 
velocity  and  pressure  fields  agree  exactly  with  known 
analytical  expressions  and  the  question  to  be  answered 
concerns  the  behavior  of  the  fiber  orientation  response. 
Figure  5.3  displays  the  angle  of  orientation  calculated 
along  the  same  path  line  per  each  respective  grid.  The 
fibers  were  initially  aligned  perpendicularly  to  the  flow 
direction  and  the  orientations  shown  in  the  table  were 
assumed  as  a  function  of  the  axial  coordinate.1 

JThe  process  integrates  along  path  lines,  but  for  the 
case  of  Poiseuille  flow  the  path  lines  are  parallel  to 
the  axial  coordinate. 


Figure  5.2  Respective  grids  for  conducting  a 
convergence  study 
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Next,  it  is  desirable  to  approximate  the  order 
of  convergence  from  the  data  of  Figure  5.3.  Initially, 
suppose  that  the  truncation  error  at  each  step  along  the 
streamline  can  be  approximated  by  the  form 

^  ^  exact  '  ♦<*> computed  “  C<h)P  <5' 

where  C  and  p  are  constants  and  p  is  the  order  of  the 
method.  Also,  h  is  defined  to  be  the  step  size  (depen¬ 
dent  on  the  mesh  refinement) ,  and  the  left-hand  side  of 
eq.(5.1)  is  the  discretized  error. 

Consider  the  data  labeled  "i"  in  going  from  the 

meshes  of  Figure  5.2a  to  5.2b. 

-0.599884  +  0.616629  -  C  (  l/3 )  P 
-0.599884  +  0.600605  a  C(y5)P 

Dividing  these  results  and  solving  for  p  yields 

23.2  *  (5/3)P 
or  p  =  6 . 1 

For  comparison,  the  data  labeled  "k"  yields 

-0.939547  +  0.954663  *  C(y3)P 
-0.939547  +  0.940057  -  C(ys)p 

Thus, 

29.6  a  (S/3)P 
or  p  =  6.6 

and  one  concludes  that  the  method  is  of  at  least  sixth 
order  which  offers  swift  convergence.  A  brief  glance 
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Figure  5.3  Results  of  the  convergence  study 
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B 
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NOTE:  All  numbers  indicate  values  of  orientation  along  the 

same  streamline,  tp  =  -0.8  . 
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at  the  data  labeled  "in"  will  show  that  a  further  mesh 
refinement  from  Figures  5.2b  to  5.2c  is  not  needed  since 
it  appears  that  we  are  incurring  errors  due  to  the  pre¬ 
cision  of  the  method.  Only  six  significant  digits  were 
carried  throughout  the  calculations  and  the  inaccuracies 
are  most  probably  due  to  round-off  error.  If  disturbed 
by  this  error,  it  is  encouraging  to  note  the  error  in  the 
solution  of  grid  B  is  less  than  a  fraction  of  a  single 
percent,  and  for  the  current  work  is  more  than  sufficient. 

Thus,  we  have  shown  the  method  of  calculating 
fiber  orientation  to  be  rapidly  convergent  in  a  manner 

g 

where  the  discretization  error  diminishes  as  0(h)  as 
h  -*■  0  .  It  is  thus  established  that  the  method  will  con¬ 
verge  with  increased  discretization.  Further  investiga¬ 
tion  has  clearly  dictated  the  fact  that  a  highly  refined 
mesh  is  not  always  the  optimum  situation.  A  definite 
trade-off  exists  between  increased  accuracy  and  additional 
computation.  Reasonable  engineering  discretion  and  exper¬ 
ience  are  needed  to  devise  an  optimum  mesh  for  computing 
fiber  orientation. 

5 . 3  Defining  an  Orientation  Parameter 

In  an  effort  to  simulate  the  actual  flow  pattern 
observed  in  real  molding  operations,  it  was  necessary  to 


44 


follow  the  orientations  of  many  particles.  Suppose  one 
designates  10  streamlines  of  a  flow  field  along  which 
one  wishes  to  monitor  fiber  orientation.  On  the  average, 
the  orientation  shall  be  recorded  approximately  20  times 
per  streamline.  Recall  that  the  value  for  the  numeri¬ 
cally  integrated  fiber  orientation  is  found  at  Pq  ,  P^  , 

•  •  •  '  Pjg  '  where  N  equals  the  number  of  intersections 
between  streamlines  and  element  boundaries.  In  addition, 
if  we  propose  5  initial  orientations,  then  the  total  num¬ 
ber  of  discrete  fiber  orientations  is  1,000  .  That  is  to 
say,  there  will  be  1,000  little  fibers  plotted  at  their 
respective  orientation  and  position  within  the  flow  field. 
While  qualitatively  appealing,  this  procedure  for  repre¬ 
senting  the  fiber  orientations  becomes  cumbersome.  Ulti¬ 
mately  one  seeks  a  more  feasible  way  of  interpreting  the 
fiber  orientation  data. 

It  is  desired  to  follow  a  development  that  appeals 
to  the  experimentalist.  Consider  a  region  of  the  flow 
field  in  which  the  fibers  are  oriented  arbitrarily  as 
shown  in  Figure  5.4.  Initially  one  defines  a  set  of  ini¬ 
tial  coordinate  axes  (x,y) .  With  respect  to  this  coordinate 
system,  the  arithmetic  mean  value  of  the  various  orienta¬ 
tions  is  calculated  and  denoted  as  lj>.  Individual  fiber 
orientations  with  respect  to  !j>  are  represented  by  and 
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Figure  5.4  Determining  the  local  fiber 

orientations  and  their  relation 
with  the  global  coordinate  system 


are  shown  in  Figure  5.4.  From  these  quantities  one  may 
define  a  root  mean  square  (rms)  value  of  the  orientation 
angle  given  by  eq.(5.2) 


$ 

rms 


(5.2) 


Clearly,  due  to  sumraetry,  the  permitted  values  for  $  & 

reside  between  0  and  ir/2.  The  symmetry  argument  stems 
from  the  fact  that  <{k  =  ($  +  tt/2)^  and  -<Jk  =  .  To  further 

simplify  the  orientation  idea,  one  normalizes  the  value 
of  $rms  by  introducing  the  orientation  parameter  f,  defined 
by  eq.  C5. 3) 
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2 

f  =  2  (cos  <j>  )  -  1  (5.3) 

Yrms 

The  values  of  f  range  from  -1  to  1  and  have  the  following 
connotations : 

-1  =>  perfect  alignment  of  fibers  transverse  to  d> 

0  =>  random  alignment 

1  =>  perfect  alignment  of  fibers  parallel  to  !j) 

Now,  having  the  interpretation  of  the  orientation  parameter 
(f)  in  hand,  one  may  graphically  display  the  orientation 
surface,  f(x,y),  within  the  domain.  Also,  quite  often  con¬ 
tours  of  f  are  plotted  to  describe  the  relative  degree  of 
fiber  orientation  that  exists  throughout  the  region. 


6.  SELECTED  NUMERICAL  EXAMPLES 


6.1  2-D  Poiseuille  Flow 


In  section  A2.8,  the  numerical  solution  for  two- 
dimensional,  pressure  driven,  shearing  flow  is  presented. 
Here,  the  results  from  the  program  COLINS  are  presented 
to  show  the  close  agreement  with  the  analytically  predicted 
orientations.  Figure  6.1a  shows  the  resulting  fiber  orien¬ 
tations  for  this  flow.  Four  initial  fiber  orientations  at 
each  point  at  the  left-hand  side  were  prescribed  to  be  0, 
±45,  90,  in  an  effort  to  simulate  a  random  distribution 
of  fibers  at  the  inlet.  This  is  to  say  that,  program  COLINS 
was  executed  four  times  with  one  of  the  specified  initial 
orientations  per  each  execution.  The  lucid  indication  of 
the  orientation  mechanisms  is  borne  out  in  this  figure.  In 
regions  near  the  exterior  boundary,  the  fibers  are  aligned 
parallel  to  the  streamlines  (see  Figure  A2.20e)  regardless 
of  their  initial  inlet  orientations .  These  were  the  same 
findings  of  the  analytic  study  of  section  4.2.1.  Contours 
of  the  orientation  parameter  are  plotted  in  Figure  6.1b, 
which  were  calculated  according  to  eq.(5.3).  The  values 
of  this  parameter  are  calculated  by  COLINS  subsequent  to 
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(a)  Numerically  predicted  fiber  orientations 
resulting  from  wall-bounded  shearing  flow 
(Four  initial  orientations  (random  dis¬ 
tribution)  at  inlet) 


Figure  6.1  Fiber  orientation  in  two-dimensional 
Poiseuille  flow 
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the  solution  of  the  individual  fiber  distributions. 
Comparing  Figure  6.1b  with  the  previously  calculated  ana¬ 
lytical  contours  of  D  (Figure  4.4a) ,  one  finds  a  close 
agreement  in  the  behavior  of  the  orientation  for  this  wall- 
bounded  shearing  flow.  Note  that  equivalence  of  the  para¬ 
meters  D  and  f  is  not  claimed,  but  merely  the  fact  that 
both  give  rise  to  the  boundary  layer  concept.  In  each  case, 
a  boundary  layer  of  oriented  fibers  near  the  wall  is  pre¬ 
dicted,  and  the  boundary  layer  thickness  is  seen  to  asymp¬ 
totically  approach  the  channel  width  at  large  values  of 
axial  coordinate. 

The  purpose  of  choosing  this  example,  once  again, 
was  to  provide  consistency  among  results,  as  the  orienta¬ 
tion  schemes  have  been  further  extended.  By  maintaining 
a  constant  check  upon  the  extensions  to  the  model,  a  high 
level  confidence  is  achieved  in  the  model  itself.  In  the 
next  sections,  another  example  to  further  test  the  numer¬ 
ical  model's  ability  to  accurately  predict  the  resulting 
fiber  orientations  in  a  shear  thinning  fluid  is  examined. 
Finally,  an  example  of  expansion  flow  and  the  resulting 
fiber  orientations  is  explored.  Here  no  analytical  solu¬ 
tion  exists  for  the  fiber  orientations  and  the  goal  of 
modeling  fiber  orientations,  numerically,  is  achieved. 
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6 . 2  2-D  Power  Law  Flow 

6.2.1  Solution  of  the  Problem 

The  problem  of  a  power  law  fluid  flowing  between 
two  parallel  walls  becomes  interesting  in  itself  when  a 
numerical  solution  is  sought.  The  finite  element  scheme 
used  previously  may  be  employed  again  with  one  major  mod¬ 
ification.  The  viscosity  of  the  fluid  is  a  function  of 
the  deformation  rate  tensor,  hence  its  value  is  not  a 
constant.  This  suggests  an  iterative  type  solution  of 
the  problem  with  one  possible  approach  outlined  below. 

Initially  the  Newtonian  velocity  field  is  obtained 
from  the  finite-element  method,  assuming  y  is  constant. 
Subsequent  to  the  determination  of  the  initial  velocity 
field,  the  viscosity  variation  is  calculated  in  terms  of 
the  spacial  coordinates.  This  is  accomplished  by  computing 
the  components  of  the  rate  of  deformation  tensor  from  the 
calculated  Newtonian  velocity  field.  The  solution  for  the 
velocity  field  is  then  repeated  using  the  variable  viscos¬ 
ity  function,  and  as  before  the  viscosity  function  is  up¬ 
dated  by  recomputing  the  components  of  the  deformation  rate 
tensor  in  terms  of  the  iterated  velocity  field.  This  pro¬ 
cedure  stops  when  the  numerical  values  for  the  components 
of  the  velocity  field  are  sufficiently  close  to  the  previous 


ones.  Now  assuming  that  the  method  of  solution  is  well 
in  hand,  one  turns  to  the  task  of  choosing  appropriate 
boundary  conditions. 

6.2. 1.1  Boundary  Conditions 

Figure  6.2a  illustrates  the  appropriate  mesh  pat¬ 
tern  for  the  two-dimensional  power  law  flow,  taking  advan¬ 
tage  of  the  existing  symmetry.  The  boundary  conditions 
along  the  solid  wall  and  channel  centerline  are  trivially 
applied  in  each  case.  The  "no  slip"  condition  at  the  wall 
is  satisfied  by  specifying  zero  velocity  there.  Since  the 
resulting  velocity  profile  is  assumed  to  be  symmetric  about 
the  channel  centerline,  the  shear  component  of  the  traction 
vector  vanishes  on  the  centerline.  In  addition,  the  trans¬ 
verse  velocity  (v)  is  prescribed  to  be  zero  everywhere  on 
the  boundary.  There  is  yet  to  specify  two  more  boundary 
conditions  to  the  domain:  one  each  at  the  upstream  and 
downstream  ends  of  the  region.  At  these  boundaries,  one 
may  either  specify  the  axial  component  of  the  velocity  or 
the  normal  component  of  the  traction  vector.  For  this  ex¬ 
ample,  it  has  been  chosen  to  prescribe  the  normal  component 
of  the  traction  vector,  F  .  Owing  to  the  fact  that  the 

normal  deviatoric  stresses  are  zero  for  this  flow,  F  is 

x 

identical  to  the  pressure.  One  sets  F  =0  at  the  downstream 

x 

boundary  and  hence  only  needs  to  calculate  the  pressure 
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difference,  which  shall  be  applied  on  the  boundary  at 
x  =  0  . 


Since  F  being  prescribed  actually  entails  a  prescribed 

X 

pressure  loading,  it  is  necessary  to  determine  the  nodal 
contributions  from  the  distributed  pressure  loading.  Re¬ 
call  that  the  nodal  forces  are  given  by 


j1  =  -/W.t  3S 

X^  X  X 


(6.1) 


where  Fx^  =  i-th  nodal  force, 

=  i-th  shape  function, 
t  =  surface  loading. 

This  equation  may  be  recast  graphically  as  shown  in 
Figure  6.3.  Carrying  out  the  integration  of  eq.(6.1)  and 


Figure  6.3  Representation  of  the  pressure 

distribution  by  concentrated  nodal 
forces  for  quadratic  interpolation 
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assuming  quadratic  interpolation  functions,  the  concen¬ 
trated  nodal  forces  are  seen  to  vary  as  shown  in  Figure 
6.3.  With  the  calculation  of  the  upstream  boundary  con¬ 
dition,  the  boundary  value  problem  becomes  completely  and 
uniquely  specified.  Numerical  results  from  the  solution 
of  the  problem  are  given  in  Figures  6 . 3b  -  d  and  may  be 
contrasted  with  those  of  2-D  Poiseuille  flow  in  Figures 
A2.20b-  f.  One  distinct  difference  between  the  flows  is 
seen  in  the  nature  of  the  vorticity  function,  which  is 
seen  to  be  confined  to  a  narrow  region  along  the  wall  for 
the  power  law  fluid.  This  is  characteristic  of  the  shear 
thinning  flows  and  can  be  explained  by  referencing  the 
pseudoplastic  velocity  profiles  of  Figure  A1.12. 

6.2. 1.2  Numerically  Predicted  Fiber  Orientations 

Here,  again  one  notes  the  fiber  orientations  pre¬ 
dicted  by  the  model  and  shown  in  Figure  6.4a.  The  power 
law  index  was  taken  to  be  n  =  .  2  and  a  random  orientation 
of  fibers  was  prescribed  at  the  left-hand  side  (inlet) . 
The  figure  shows  the  individual  fibers  plotted  within  the 
flow  region  and  it  is  interesting  to  note  the  expanded 
core  region  possessing  random  orientation  even  at  dis¬ 
tances  away  from  the  inlet.  Qualitatively,  this  behavior 
could  be  anticipated  from  the  information  describing  the 


flow  field. 


An  equally  interesting  observation  is  provided  by 
Figure  6.4b,  where  the  reduction  in  the  boundary  layer  is 
associated  with  increasing  shear  thinning.  The  boundary 
layer  for  this  example  has  been  taken  as  f  j>  .  8  which 
corresponds  to  a  highly  oriented  region  of  fibers  paral¬ 
lel  to  the  wall.  As  an  added  note,  it  can  be  shown  ana¬ 
lytically  that  the  boundary  layer  will  approach  the  chan¬ 
nel  half— width  as  x-*-°°.  The  degree  of  pseudoplasticity 
will  determine  the  order  of  the  singularity  at  y  =  a  and 
hence  how  fast  the  boundary  layer  approaches  the  channel 
half-width.  This  point  was  also  made  mention  of  in 
section  4.2.2. 

6 . 3  Expansion  Flow 

Having  thus  far  gained  considerable  confidence 
in  the  numerical  scheme  for  predicting  fiber  orientations, 
one  arrives  at  a  flow  configuration  for  which  nonanalyt- 
ical  solutions  exist.  It  becomes  necessary  to  solve  the 
problem  entirely  by  numerical  means.  This,  of  course, 
has  been  our  ultimate  objective  all  along,  but  it  was 
important  to  check  the  method  against  known  solutions  at 
each  stage  in  the  development  of  the  model. 

Flow  through  a  10:1  expansion  is  investigated 
in  this  section  to  view  the  fiber  orientations  in  an 


(a 


)  Numerically  predicted  fiber  orientations 
in  a  steady  power  law  flow  (n=.2) 

(Four  initial  orientations  (random  dis¬ 
tribution)  at  inlet) 


(b)  Contours  of  orientation  parameter  (f=.8) 
for  varying  power  law  index 


Figure  6.4  Fiber  orientations  in  two-dimensional 
power  law  flow 
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extensional  flow  situation.  While  the  flow  itself  is 
not  purely  elongational  the  region  removed  from  the 
immediate  entrance  consists  mainly  of  an  extensional 
flow.  For  this  reason  the  dominant  behavior  of  the  nor¬ 
mal  velocity  gradients  may  not  be  ignored.  The  flow 
configuration  under  present  consideration  is  encountered 
in  practice  at  the  entrance  to  a  large  cavity.  As  usual 
we  are  assuming  sufficiently  slow  motions  so  that  the 
Stoke' s  equations  of  motion  are  valid.  This  in  essence 
assumes  that  all  transient  and  inertia  dominated  behavior 
is  neglected.  The  fluid  is  considered  to  behave  in  the 
Newtonian  sense,  with  the  extension  to  power  law  fluids 
following  immediately  in  light  of  the  previous  discussion 
of  section  6.2. 

6.3.1  Solution  of  the  Problem 

6.3. 1.1  Boundary  Conditions 

The  appropriate  boundary  conditions  for  this  prob¬ 
lem  are  labeled  in  Figure  6.5,  which  also  illustrates  the 
grid  design  for  the  problem  solution.  The  boundary  con¬ 
ditions  by  now  are  rather  obvious  with  the  possible  excep¬ 
tion  perhaps  of  the  far  field  boundary  condition  applied 
at  r=10.  Again,  as  before,  one  has  the  freedom  of  speci¬ 
fying  force  components  or  velocity  components  on  this 


59 


portion  of  the  domain.  It  is  desired  to  prescribe  the 
velocity  components  at  far  field  since  they  may  be  cal¬ 
culated  exactly  from  the  solution  to  the  problem  of 
Jef fery-Hamel  flow.  For  source  flow  in  a  diverging 
section  at  small  Reynold's  numbers  the  velocity  profile 
for  Jef fery-Hamel  flow  has  been  shown  to  be  [35] 


•  2„ 

u  _  1  sin  9 

u  .2 

o  sm  a 


(6.2) 


where  u  =  4/rtr  for  unit  flow  rate  and  a  =  ir/2.  The  u 
o  o 

is  determined  from  continuity  considerations  in  much  the 
same  manner  as  v  is  determined  in  eq.(A1.64) .  Note 

aVy 

that  this  velocity  profile  is  parabolic  in  the  0-direction 
and  also  maintains  a  radial  dependence  through  the  constant 


Figure  6.6  Jef fery-Hamel  flow  in  a 
diverging  section 
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.  At  large  radius  values  for  our  expansion  flow  problem, 
one  may  specify  the  Jef fery-Hamel  velocity  profile  (eq. 

(6.2))  with  good  accuracy.  It  is  assumed  that  when  r =  10 
the  effects  of  the  entrance  region  are  negligible  and  the 
flow  at  this  point  is  as  though  it  emanated  from  a  source. 
This  was  the  prime  consideration  for  choosing  a  large 
value  of  r  at  which  to  apply  the  final  boundary  condition. 

The  solution  to  the  fluid  mechanics  problem  is 
displayed  in  Figures  6 . 7a  -  d  in  the  usual  spirit  of  con¬ 
tour  plots.  Unless  otherwise  noted,  ten  contours  are  plot¬ 
ted  per  each  field  variable  equally  spaced  between  the  max¬ 
imum  and  minimum  values.  Contours  that  are  not  particularly 
smooth  are  a  result  of  plotting  a  field  variable  that  has 
a  value  near  zero.  Clearly,  if  the  noise  level  of  the  solu¬ 
tion  is  even  minimal,  it  will  be  reflected  in  such  a  con¬ 
tour.  In  any  event,  perusing  Figures  6.7a-  d,  one  may  be 
convinced  that  the  solution  to  the  fluid  mechanics  problem 
is  a  correct  one. 

6.3. 1.2  Numerically  Predicted  Fiber  Orientations 

In  all  the  previous  examples  and  discussion,  it 
has  been  shown  how  fibers  align  parallel  to  streamlines 
under  steady  shearing  flow.  So  far,  in  the  present  dis¬ 
cussion  one  has  dealt  with  the  flow  in  an  expansion  and 


MAX  =  1.07  MIN  =  -4.06 


Figure  6.7a  Shear  stress  contours  of  the  expansion 
flow  problem 


Figure  6.7b  Contours  of  axial  veloctiy  of  the  expan¬ 
sion  flow  problem 


MAX  =  0.0  MIN  =  -1.0 


Figure  6.7c  Contours  of  the  stream  function  of  the 
expansion  flow  problem 


Figure  6.7d  Contours  of  vorticity  of  the  expansion 
flow  problem 
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now  inquires  as  to  the  resulting  fiber  orientations 
from  such  a  flow.  Obviously,  from  the  knowledge  of  the 
qualitative  orientation  mechanisms  from  section  3.1, 
one  would  expect  to  find  transverse  orientations  of  the 
fibers  in  the  region  of  the  flow  where  elongational  be¬ 
havior  dominates.  Figure  6.8  shows  exactly  this  phenom¬ 
enon.  Fibers  are  completely  aligned  in  the  constricted 
region  owing  to  the  high  vorticity  and  shearing  compo¬ 
nents  existing  there.  In  practice  this  region  may  corres¬ 
pond  closely  to  the  behavior  experienced  in  the  sprue 
of  the  mold.  As  the  fibers  enter  the  mold  cavity  they 
are  oriented  transversely  to  the  streamlines  primarily 
due  to  the  elongational  component  of  the  flow.  Another 
interesting  point  is  to  note  the  shearing  mechanism  for 
orientation  at  work  all  along  the  wall. 

Figure  6.9  displays  the  orientation  parameter  for 
the  expansion  flow  problem.  Recall  that  a  contour  value 
of  1.0  corresponds  to  a  completely  aligned  configuration. 
The  angle  of  alignment  with  the  global  axes  may  be  found 
in  determining  the  quantity  <P .  (This  quantity  was  defined 
in  section  5.3.}  Alternatively,  one  may  deduce  the  steady 
state  angle  of  alignment  from  a  knowledge  of  the  local 
flow  kinematics.  That  is  to  say,  if  the  local  flow  char¬ 
acteristics  are  predominately  shearing,  then  f  =  1.0 
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corresponds  to  aligned  fibers  parallel  to  the  stream¬ 
lines.  Transverse  orientations  to  the  fluid  steamlines 
at  steady  state  result  from  regions  of  the  flow  dominated 
by  elongational  components . 


7.  EXTENSIONS  TO  THE  THEORY 


Throughout  this  documentation,  mention  has  been 
made  of  several  shortcomings  to  the  fiber  orientation 
theory  as  presented  herein.  While  it  has  been  shown  to 
provide  excellent  results  within  the  assumptions  set 
forth,  there  are  a  few  extensions  to  the  existing  ideas 
that  need  to  be  brought  to  light. 

The  extension  of  the  model  to  include  axisymmetric 
flows  is  straightforward  in  theory.  It  should  be  realized 
that  both  the  finite-element  method  and  the  orientation 
equations  (eqs . (4 . 1-2) )  need  to  be  generalized  to  handle 
axisymmetric  flows.  Neither  is  suspected  of  causing  any 
conceptual  difficulties  since  only  a  coordinate  transfor¬ 
mation  is  encountered.  For  example,  eq.(4.3)  may  be  writ¬ 
ten  in  terms  of  axisymmetric  components  by  the  usual  rules 
of  tensor  transformations.  Hence,  the  present  analysis 
could  be  extended  to  include  axisymmetric  flows  directly. 

Extending  the  analysis  of  the  model  to  a  full  three 
dimensions  again  provides  little  difficulty  conceptually. 
The  limitation  of  this  extension  usually  becomes  apparent 
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in  the  solution  to  the  finite  element  problem.  What 
usually  happens,  using  three-dimensional  elements,  results 
in  the  rapid  exhaustion  of  computer  core  storage.  As  a 
guideline,  the  consideration  of  a  three-dimensional  domain 
increases  the  number  of  unknowns  by  a  factor  of 

2  (94/33)  *n3  -  0(6n3) 

where  n3  is  the  number  of  elements  through  the  thickness 
of  the  domain.  The  ratio  of  (94/33)  is  the  increase  in 
the  number  of  degrees  of  freedom  in  going  from  the  6  node 
triangular  element  to  the  10  node  tetrahedral  element. 

Here,  one  considers  the  velocities,  stresses  and  pressure 
as  primitive  variables.  Indeed,  there  are  methods  of  solu¬ 
tion  which  reduce  the  overall  core  storage  requirements 
(see  section  A2.6)  but  considerable  effort  would  be  ex¬ 
pended  to  develop  the  necessary  software.  In  this  spirit, 
the  extension  to  three-dimensional  domains  seems  somewhat 
less  tractable  than  before. 

It  has  been  suggested  that  many  of  the  practical 
molding  operations  poorly  approximate  steady  state  flow 
kinematics,  especially  in  the  proximity  of  the  flow  front. 
Tadmor  [36]  suggests  that  the  flow  front  actually  plays  a 
significant  role  in  the  orientation  of  fibers.  For  these 
reasons,  it  would  be  beneficial  to  include  the  transient 
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terms  of  the  equations  of  motion  in  the  solution  of  the 
flow  problem.  This  additional  modification  would  be  es¬ 
pecially  helpful  if  it  were  developed  using  a  Lagrangian 
method  of  solution.  Here,  the  grid  is  allowed  to  deform 
continuously  to  monitor  the  actual  behavior  of  the  ad¬ 
vancing  free  surface.  Thus  the  mold  filling  process  could 
be  closely  modeled  by  marching  out  the  solution  in  the 
time  domain  with  continuous  mesh  updating  and  rearrange¬ 
ment.  At  this  point,  it  is  not  known  how  well  Jeffery's 
theory  of  particle  orientations  would  approximate  the  true 
behavior  of  the  solution,  since  the  inclusion  of  the  tran¬ 
sient  terms  has  violated  Jeffery's  assumption  of  steady 
flow.  This  ramification  needs  to  be  determined. 

The  important  consideration  of  including  the  ef¬ 
fect  of  fiber  orientation  on  the  structure  of  the  velocity 
field  has  been  discussed  in  section  5.1.  This  iterative 
procedure  seems  highly  feasible  once  the  orientation- 
dependent  viscosity  is  known.  Recent  work  on  obtaining 
such  viscosity  functions  [24]  has  shown  that  this  effect 
may  be  important  for  certain  flow  situations.  Clearly, 
this  is  an  avenue  for  further  research. 

If  the  constitutive  equation  describing  the  flow 
includes  the  possibility  for  the  fluid  to  possess  some 
elasticity,  a  viscoelastic  model  must  be  employed  to 
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determine  behavior  of  the  flow.  The  numerical  solution 
of  the  viscoelastic  fluids  problem  can  be  accommodated 
by  existing  computer  codes  if  the  elasticity  effects 
within  the  fluid  are  not  too  great.  In  addition,  one 
other  modification  needs  to  be  made.  The  Jeffery  equa¬ 
tions  (eqs .  (4 . 1-2) )  that  describe  particle  orientations 
in  viscoelastic  fluids  must  be  re-evaluated  in  light  of 
the  non-Newtonian  fluid  behavior.  These  modified  equa¬ 
tions  may  be  referenced  in  Van  Oene  [39],  for  example, 
and  the  incorporation  of  this  extension  seems  straight¬ 
forward,  at  least  in  principle. 

Concluding  this  discussion,  then,  a  final  comment 
is  in  order.  Several  possible  extensions  to  the  orien¬ 
tation  theory  have  been  suggested,  whose  effects  on  the 
fiber  orientations  are  unknown.  It  is  possible  that  the 
inclusion  of  the  proposed  extensions  will  not  significantly 
alter  the  theoretical  results  in  light  of  the  excellent 
findings  predicted  by  the  current  theory.  The  empirical 
and  analytical  verification  to  the  present  theory  yields 
a  high  level  of  confidence  in  predicting  the  fiber  orien¬ 
tations  in  real  flows,  and  a  more  sophisticated  model  may 
not  deliver  additional  insight. 


CONCLUSIONS 


A  method  has  been  developed  to  determine  the 
fiber  orientations  in  a  dilute  suspension  undergoing 
arbitrary,  two-dimensional  deformations.  For  flow  situ¬ 
ations  where  shearing  deformation  dominates ,  the  dynamic 
orientation  of  the  fibers  has  been  predicted.  The  equil¬ 
ibrium  orientation  for  shearing  flows  has  been  shown  to 
be  parallel  to  fluid  streamlines.  For  cases  of  wall- 
bounded  channel  flow,  this  orientation  mechanism  is  re¬ 
sponsible  for  the  occurrence  of  a  boundary  layer  phenom¬ 
enon  that  produces  regions  of  highly  aligned  fibers  along 
the  channel  walls.  Perturbations  from  the  Newtonian  con¬ 
stitutive  assumption  are  possible  by  considering  the 
analysis  to  be  quasi-Newtonian.  Here,  the  effect  of  in¬ 
creasing  pseudoplasticity  reduces  the  thickness  of  the 
boundary  layer. 

For  certain  problems  where  analytical  expressions 
do  not  exist  for  the  field  variables,  the  finite-element 
method  is  employed.  Upon  solution  to  the  fluid  mechanics, 
the  orientation  equation  may  be  numerically  integrated 
along  a  path  line  to  yield  the  resulting  rotations.  This 
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methodology  has  been  incorporated  into  a  computer  routine 
that  allows  one  to  calculate  the  resulting  fiber  orien¬ 
tations  for  rather  complex  geometry  flows. 

Employing  the  above  methodology  to  the  problems 
of  Newtonian  flow  in  the  vicinity  of  an  infinite  expan¬ 
sion,  the  resulting  equilibrium  orientations  have  been 
shown  to  be  transverse  to  the  fluid  path  lines.  These 
general  trends  in  fiber  orientations  are  also  expressed 
intuitively  by  considering  the  stresses  acting  on  a  sin¬ 
gle  fiber  immersed  in  a  flow  of  given  kinematics. 
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APPENDIX  1 


FLUID  MECHANICS 


In  this  chapter  the  basic  conservation  concepts 
that  provide  a  foundation  to  the  fluid  mechanics  problem 
are  examined.  The  equations  that  are  discussed  are  gener¬ 
ic  in  nature,  since  they  apply  to  all  material  systems. 

The  question  of  the  behavior  of  various  fluids  falls  under 
the  category  of  constitutive  relationships  to  be  discussed 
in  section  A1.2.  It  is  the  coupling  of  both  conservation 
equations  and  constitutive  approximations  that  generate 
the  method  of  solution  to  practical  engineering  problems . 

A  Cartesian  reference  frame  has  been  adopted  to 
develop  the  governing  equations.  For  conciseness,  indicial 
or  Einsteinian  notation  will  be  introduced  throughout  this 
discussion,  with  the  usual  assumption  that  repeated  indices 
are  summed  over  the  range  1-3,  unless  otherwise  noted.  A 
more  elegant  derivation  utilizes  general  tensorial  mathe¬ 
matics  to  develop  a  set  of  governing  equations  insensitive 
to  any  particular  coordinate  system.  The  application  to  a 
useful  reference  frame  is  then  only  an  exercise  in  manipula¬ 
tion  to  convert  tensorial  quantities  to  physical  quantities. 
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This  line  of  approach  tends  to  be  somewhat  abstract  and 
physical  insight  is  often  overshadowed. 

Often,  in  an  effort  to  grasp  new  ideas,  it  is 
instructive  to  relate  them  with  a  more  familiar  area. 
Recall,  from  linear  elasticity,  that  the  general  formu¬ 
lation  of  a  problem  requires  15  independent  equations  to 
completely  specify  the  15  unknowns.  The  unknowns  are  6 
stress  components,  6  strain  components,  and  3  displacement 
field  components ;  and  the  equations  relating  these  unknowns 
are  3  equilibrium  equations,  6  constitutive  equations,  and 
6  compatibility  conditions.  Thus,  neglecting  the  details 
of  the  boundary  conditions,  this  set  of  15  equations  and 
15  unknowns  forms  the  basis  of  a  well-posed  problem  for 
linear,  isothermal  elasticity. 

Alternatively,  one  shall  see  that  in  the  fluid 
mechanics  formulation,  that  the  total  number  of  unknowns 
is  10:  6  stress  components,  3  velocity  components,  and 

the  density.  In  order  for  this  problem  to  be  well-posed 
it  is  essential  that  10  independent  equations  exist.  Con¬ 
servation  of  mass  provides  one  equation,  while  the  con¬ 
servation  of  linear  momentum  adds  three  more  independent 

relations.  So,  it  is  apparent  that  to  solve  fluid  mechan¬ 
ics  problems,  one  must  introduce  six  additional  relation¬ 
ships.  Some  assumptions  may  result  from  physical  insight 
into  the  problem  that  allows  an  unknown  to  be  set  equal  to 
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zero.  Constitutive  approximations  introducing  relation¬ 
ships  among  stresses  and  strain  rates  also  allow  for  addi¬ 
tional  independent  equations.  As  previously  noted,  these 
will  be  investigated  in  section  A1 . 2 . 

A1 . 1  General  Concepts 

A1 . 1 . 1  Conservation  of  Mass 

As  for  any  continuum,  a  volume  of  fluid  must  sat¬ 
isfy  the  basic  law  of  conservation  of  mass,  which  simply 
states  that  all  mass  within  a  given  control  volume  must  be 
accounted  for  at  all  times.  Added  to  this  is  the  assump¬ 
tion  of  no  atomic  disintegration  taking  place  within  the 
specified  control  volume.  Consider  the  control  volume 


m2 

Figure  A1 . 1  Control  volume  for  conservation 
of  mass 


81 


Applying  the  conservation  of  mass  principle,  one  can  write 
in  a  statemental  form 

rate  of  mass  -  rate  of  mass  =  rate  of  mass  (Al.l) 

into  c.v.  out  of  c.v.  accumulation 

within  c.v. 

A  suitable  expression  for  the  rate  of  mass  (m)  will  be 
given  in  general  by  the  product  of  the  density  and  the 
flow  rate,  as  given  by  eq.(A1.2), 

m  =  pvdA  ( A1 . 2) 

where  p  is  defined  to  be  the  fluid  density  and  (vdA)  is 
the  corresponding  flow  rate  through  a  given  face  of  the 
control  volume.  Utilizing  the  left-hand  side  of  eq.(Al.l) 
in  the  x^  direction  yields 

m^  -  (m1  +  dx1)  ,  (A1.3) 

and  per  eq. (Al. 2) 

- (pv1dx2dx3)  dx1  =  (pv1)dV  (A1.4) 

If  one  uses  the  same  procedure  for  the  remaining  coordinate 
directions,  it  is  apparent  eq.(Al.l)  becomes 

(pv3)JdV  =  (|^)dV 

(A1.5) 

and  since  the  control  volume  has  been  chosen  arbitrarily 


(pv1)  - 


(Pv2)  - 
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eq.(A1.5)  can  be  simplified  to 

9  p  ,  3  ,  . 

at  +  37"  (pVi}  =  0  ( A1 . 6 ) 

l 

At  once,  one  can  see  the  compactness  of  eq.(A1.6)  in  its 
indicial  form.  Recall  from  the  previous  conventions  that 
the  second  term  of  eq.(Al.6)  implies  three  individual 
terms  (the  index  i  is  summed  from  1-3).  Eq.(A1.6)  is 

the  most  general  expression  for  the  conservation  of  mass. 
Some  textbooks  also  refer  to  this  equation  as  the  conti¬ 
nuity  equation  in  that  it  assures  compatibility  within 
the  flow  domain.  Either  expression  is  correct  and  both 
shall  be  utilized  throughout  this  discussion. 


Looking  at  eq. (A1.6) ,  it  may  be  possible  to  make 
some  simplifying  assumptions.  If  the  analysis  is  one  of 
steady  state,  this  implies  that  the  time  dependent  term 
disappears  from  eq.(A1.6)  and  one  is  left  with 


(pv^ 


0  . 


(A1.7) 


Next,  if  it  is  assumed  that  the  density  variation  with 
respect  to  coordinate  positions  is  negligible,  then  eq. 
(A1.7)  reduces  to 


3x . 

l 


0 


(Al.  8) 
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The  last  assumption  was  one  of  incompressibility,  a  be¬ 
havior  to  which  many  engineering  fluids  closely  approx¬ 
imate.  Again  one  should  be  cautioned  that  eq.(A1.8)  ap¬ 
plies  only  to  Cartesian  coordinate  systems,  but  trans¬ 
formations  to  other  reference  frames  are  found  in  the 
literature. 

Al.1.2  Conservation  of  Linear  Momentum 

If  one  applies  Newton's  second  law  (eq.(Al.9)  in 
inertial  coordinates  to  the  same  control  volume  as  before, 
shown  for  convenience  in  Figure  A1.2,  he  will  arrive  at 
the  expression  for  the  conservation  of  linear  momentum. 


Figure  A1.2  Control  volume  for  conservation 
of  linear  momentum 
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The  first  question  to  be  answered  concerns  the  choice  of 
coordinate  systems.  It  was  convenient  in  section  Al.1.1 
to  choose  an  inertial  frame  of  reference  to  discuss  the 
conservation  of  mass  condition.  Now  it  will  be  necessary 
to  speak  of  a  coordinate  system  imbedded  within  the  con¬ 
trol  volume  and  moving  with  the  local  velocity  of  the 
fluid.  Why  does  one  seek  a  desired  change  of  reference 
frame?  It  is  simply  for  convenience  in  describing  the 
local  accelerations  and  one  shall  see  that  this  approach 
leads  to  two  distinct  acceleration  terms . 

Since  one  has  elected  to  utilize  a  material  coor¬ 
dinate  system,  the  velocities  are  functions  of  both  space 
coordinates  and  time.  Implementing  the  use  of  indicial 
notation,  the  resultant  velocity  of  the  control  volume 
becomes 

(vi(xi,t)  +  dvi(xi,t))  -  (vi(xi,t)) 
and  per  chain  rule 

3v.  3v. 

dvitxi't)  IxJ  dxj  +  TT  dt  • 

Substituting  eq.(Al.ll)  into  eq.(A1.9),  which  is  valid 
for  inertial  systems,  yields 


dx  . 

3v.  ] 

3v.  3v.l 

k 

-*?■  + 
dt 

1 

St 

J 

dV  = 

p 

1  ,  l 

3 1  j  3x  . 

^  3  , 

=  dvi(xi,t) 

(A1.10) 


(Al.ll) 


F  . 
x 


P 


dV  .  (A1.12) 
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The  first  term  on  the  right-hand  side  of  eq.(A1.12)  is 
the  rate  of  change  in  velocity  with  respect  to  the  mate¬ 
rial  coordinate  system,  whereas  the  second  term  accounts 
for  any  accelerations  of  the  material  coordinate  system 
with  respect  to  an  inertial  reference  frame.  The  latter 
effect  is  sometimes  termed  the  transport  or  convective 
part  of  the  acceleration  and  together  with  the  local  rate 
of  change  in  velocity  represents  the  total  or  material 
derivative  of  velocity  with  respect  to  some  inertial  ref¬ 
erence  frame. 

To  clarify  the  above  discussion,  an  example  is  in 
order.  Consider  the  case  of  liquid  in  a  long  rotating 
cylinder,  where  the  external  wall  spins  with  constant 
angular  velocity.  The  most  natural  reference  frame  to 
describe  the  fluid  kinematics  is  one  of  cylindrical  coor¬ 
dinates  (r,0,z),  shown  in  Figure  A1.3. 


Figure  A1.3  Fluid  in  a  long,  rotating  cylinder 
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Since  the  flow  is  steady  by  assumption,  the  local  accel¬ 
eration  term  in  eq. (A1.12)  is  zero.  The  convective  terms, 
however,  are  nonzero  in  that  the  material  coordinate  sys¬ 
tem  is  rotating  with  respect  to  some  fixed  reference 
frame . 


Next,  one  may  apply  a  force  balance  to  the  control 
volume  of  Figure  A1 . 4  in  order  to  compute  the  forces 
acting  on  a  fluid  element.  From  the  figure  it  is  clear 
that  both  body  forces  and  surface  tractions  act  upon  the 
fluid  element.  The  body  force  may  consist  of  three  com¬ 
ponents,  each  of  which  act  along  the  respective  coordinate 
directions.  However,  the  convention  pertaining  to  surface 


Figure  A1.4  Control  volume  for  force  balance  on 
a  fluid  element  in  the  xn -direction 
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tractions  needs  to  be  elaborated  upon.  Surface  tractions 
are  defined  in  terms  of  stresses,  whose  components  trans¬ 
form  like  those  of  a  second  rank  tensor.  Apparently  two 
subscripts  are  needed  to  describe  a  given  component  of 
stress.  We  adopt  the  notation  such  that  the  first  sub¬ 
script  defines  the  surface  upon  which  the  stress  acts  (the 
surface  is  given  by  the  unit  outward  normal  vector)  while 
the  second  subscript  dictates  the  direction  of  the  stress 
component.  A  quick  glance  at  Figure  A1.4  will  illustrate 
the  use  of  this  convention. 


Force  balances  in  the  x 2  and  x^  coordinate  directions 
yield  similar  expressions  and  using  indicial  notation 
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may  be  combined  to  read 


■  3a  .  .  'i 

--J—  +  pf.  dV  =  F. 
OXj  lj  i 


Substituting  eq.(A1.14)  into  eq.(A1.12] 


(A1 .14) 


3v .  3v. 

„  r  ,  l 

P  -zrr~  +  V .  -r - 

3t  J  3x . 

^  3 


3a  .  . 

+  Pf. 

3x .  x 

3 


(A1.15) 


one  arrives  at  the  final  expression  for  the  principle  of 
conservation  of  linear  momentum. 

Al . 1 . 3  The  Stress  Tensor 

One  of  the  most  difficult  points  to  understand 
when  first  learning  fluid  mechanics  is  the  relation  be¬ 
tween  normal  stresses  and  pressure.  The  ideas  to  be  in¬ 
troduced  here  are  compatible  with  those  learned  from  an 
undergraduate  fluid  mechanics  course,  say,  and  some  com¬ 
parison  may  be  helpful. 


Thus  far,  the  stresses  that  have  been  referred 
to  were  total  stresses  (cr  )  •  The  components  of  the 
total  stress  tensor  represent  a  fully  populated  second 
rank  tensor,  shown  below.  Diagonal  components  represent 


1 
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stresses  which  tend  to  dilate  the  fluid,  while  off- 
diagonal  stresses  correspond  to  shearing  stresses. 

Recall  now  that  the  term  "pressure"  refers  to  only  dila- 
tational  behavior  and  must  therefore  be  related  uniquely 
to  the  normal  stresses  (a  ,  c?22  ,  a^)  .  Hence,  it 
seems  reasonable  to  define  the  pressure  as  the  average 
value  of  the  normal  stresses  given  by  eq.(A1.16) 

P  =  -  1/3(t11  +  a ^  +  ^33)  ( Al .  16a) 

=  -  V3  oLi  .  (Al  .16b) 

The  negative  sign  indicates  that  pressure  is  inherently 
compressive  in  nature.  Note  that  this  definition  is 
compatible  with  the  notion  of  hydrostatic  pressure  if 
the  fluid  is  at  rest.  In  this  case,  all  three  normal 
stresses  are  identical  and  their  average  defines  the 
hydrostatic  pressure. 

Next,  one  may  introduce  the  important  concept  of 
deviatoric  stress.  As  a  consequence  of  eq.(A1.16),  the 
total  stress  tensor  may  be  represented  by  two  terms: 


(Al  .17) 


(Al .18) 
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in  which  the  6^  is  defined  to  be  the  Kronecker  delta. 

It  takes  on  the  value  of  unity  for  i=j  and  is  zero  other¬ 
wise.  The  tensor  represents  that  part  of  the  stress 

tensor  as  a  result  of  fluid  motion  or  deformation.  It 
is  commonly  called  the  deviatoric  stress  tensor.  An 
alert  reader  usually  inquires  as  to  the  nature  of  the 
diagonal  terms  of  the  P^j  tensor.  Quite  simply  they  are 
the  normal  stresses  that  arise  from  deforming  the  fluid. 
They  are  not  related  to  the  pressure,  since  the  pressure 
has  been  extracted  as  a  consequence  of  eq.(A1.17).  The 
confusion  usually  stems  from  the  fact  that  a  first  course 
on  the  subject  only  emphasizes  the  shearing  behavior  of 
fluids,  say  for  instance,  water  flowing  in  a  tube.  Little 
reference  was  made  of  the  fact  that  some  fluid  systems 
exhibit  anomalous  normal  stresses.  As  an  example,  egg 
whites  hanging  from  their  shell  (Figure  A1.5)  produce 
substantial  normal  stresses  ( P 33 )  within  the  bulk  of  the 
fluid.  These  normal  stresses  are  seen  to  be  at  least 
large  enough  to  support  the  weight  of  the  fluid  in  ten¬ 
sion.  Other  such  fluids  exhibiting  this  exotic  behavior 
include  body  mucuses  and  a  slew  of  industrial  polymers. 
Even  the  familiar  concept  of  surface  tension  is  really 
another  illustration  of  normal  stresses,  that  can  develop 
under  many  circumstances. 
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Figure  A1.5  An  illustration  of  the  normal 
stresses  found  in  egg  whites 

Before  leaving  this  discussion  of  normal  stresses, 
one  last  comment  is  in  order.  When  one  experimentally  mea¬ 
sures  "pressure"  in  the  laboratory,  he  is  actually  measur¬ 
ing  a  total  stress  component  (for  definiteness  a^)  .  From 
eq.(A1.18)  it  is  apparent  that  the  total  stress  is  really 
the  sum  of  the  isotropic  pressure  and  the  deviatoric  stress 
component.  It  is  only  through  an  appropriate  rheological 
constitutive  equation  that  one  usually  claims  P  to  be 
zero.  Hence,  inadvertantly  one  is  measuring  the  quantity 
known  as  pressure,  but  the  reader  is  cautioned  that  occa¬ 
sionally  normal  deviatoric  stresses  result  from  the  fluid 
motion  and  one  no  longer  measures  just  the  "pressure". 
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Al.1.4  Conservation  of  Angular  Momentum 


Consider,  again,  the  control  volume  of  Figure  A1 . 4 . 
In  addition  to  the  linear  motion  that  these  stresses  pro¬ 
duce,  there  also  exists  angular  motion  (rotation)  that  must 
be  conserved.  The  conservation  of  angular  momentum  is  the 
analogue  of  eq.(A1.9),  with  forces  being  replaced  by  torques 
and  angular  velocities  substituted  for  linear  ones.  So, 


Figure  A1.6  shows  the  same  control  volume  as  Figure 
Al.4  with  some  of  the  stresses  deleted  for  clarity.  Also 


Figure  A1.6  Control  volume  for  moment 
balance  about  the  axis 


note  that  the  coordinate  system  has  been  placed  at  the 
center  of  the  control  volume  and  will  be  considered  as 
inertial.  In  terms  of  inertial  coordinates,  then,  con¬ 
sider  a  rotation  about  the  x^  axis  for  definiteness. 
Eq.(A1.19)  becomes  specifically 


T3  3t  {133a3] 


Substitution  for  the  moment  of  inertia  and  noting  that 
it  is  independent  of  time  yields 


t3  =  3t  (rarg3a3>  =  pdVrg3  TF 


(A1.20) 


where  rg^  is  the  radius  of  gyration  about  the  x^  axis. 

The  torque  that  the  stresses  produce  may  be  summed  as 

dx,  3a, „  dx, 

T_  =  a.  0dx_dx0 -  +  (a.,  +  -r -  dx,  )  dxodx0 - 

3  12  2  3  2  12  3x,  1  2  3  2 


-  a«,dx,dx, — : 

21  1  3  2 


So--,  dx» 

-  <°21  +  “3x7"  dx2)dxldx3— 


(A1.21) 


Note  that  the  normal  stresses  do  not  contribute  to  the 
moment  as  a  result  of  a  convenient  placing  of  the  coor¬ 
dinate  system  at  the  centroid  of  the  control  volume.  If 
one  retains  only  the  first-order  terms  of  eq.(A1.21)  and 
substitutes  into  eq.(A1.20),  the  resulting  expression  for 


conservation  of  angular  momentum  in  the  x^  direction  is 

2  ^3 

(ai2  “  a2i}  =  prg3  ~Jt  * 


(A1.22) 
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Since  eq. (A1.22)  must  hold  in  the  limit  of  the  continuum, 
the  size  of  the  control  volume  of  Figure  A1 . 6  is  allowed 


to  approach  zero.  Note  that  the  right-hand  side  of  eq. 


(A1.22)  tends  to  zero  since  r_.  ■+  0  as  c.v.  -*  0 . 

g3 

must  be  that  a12  =  a2i  for  t*1*2  equality  to  hold. 


Thus , 


it 


The  conclusion  to  be  drawn  from  this  exercise  is 
that  the  order  of  subscripts  pertaining  to  the  stresses 
does  not  matter  (i.e.  o .  .  =  ct  .  . ) .  Therefore,  the  stress 
tensor  is  seen  to  be  symmetric  in  nature  and  the  total 
number  of  unknown  independent  stress  components  reduces 
from  nine  to  six.  This  point  was  mentioned  previously; 
now  the  justification  is  established  for  such  a  claim. 


Al . 1 . 5  Material  Motion 

It  is  natural  when  referring  to  fluid  flow  to 
express  the  motion  of  the  continuum  in  terms  of  the  vel¬ 
ocity  field.  Consider  a  fluid  continuum  undergoing  an 
arbitrary  motion,  shown  in  Figure  A1.7. 


Figure  A1.7  Material  motion  for  a  fluid  continuum 
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An  inertial  Cartesian  reference  frame  is  estab¬ 
lished  to  assist  in  describing  the  particle  motions. 
Unlike  the  methodology  in  linear  elasticity,  it  is  un¬ 
desirable  to  use  displacement  gradients  in  order  to  des¬ 
cribe  material  motion.  The  primary  reason  is  that  the 
displacements  are  usually  far  greater  than  those  encoun¬ 
tered  in  solid  mechanics . 


The  difference  in  local  velocity  vectors  for  two 
neighboring  points  x0^  and  Xq^+Ax^  is  given  by  eq.(A1.23) 


vi  “  v0i 


3v^ 

dv .  =  -  dx  . 

x  ax.  i 


(A1.23) 


Again,  since  the  coordinate  system  is  fixed,  the  velocities 
are  function  only  of  position.  Clearly,  differences  in 
velocities  between  points  x0^  and  x0^+Ax^  must  introduce  a 
deformation  to  the  fluid  continuum.  Thus,  one  arrives  at 
the  conclusion  that  fluid  deformation  is  related  to  the 
velocity  gradients.  As  is  customary,  the  velocity  gradient 
tensor  may  be  written  as  two  distinct  parts  with  individual 
characteristics . 


3vi 

3v .  1 

'3vi  3v . 

y2 

9x . 

+  Zx.\ 

+  Vi 

3x .  3x. 

3 

3  1J 

(Al. 24) 


The  first  term  on  the  right-hand  side  of  eq.(A1.24)  is 
responsible  for  actual  deformation  of  a  fluid  element. 
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Hence,  it  will  be  referred  to  as  the  deformation  rate 


tensor  d^ ^ .  Upon  interchanging  i  and  j  ,  it  is  obvious 
that  this  "strain  rate"  tensor  is  symmetric. 


d.  .  =  V2 


^3v- 

3v .  • 

(  _ 

±  + 

3 

dx  . 

9x . 

V. 

3 

'll 

d12 

d13 

'12 

d22 

d23 

d  _ 

13 

23 

33 

of 

eq . 

(Al.. 

(A1.25a] 


(A1.25b; 


The  second  term  of  eq.(A1.24)  does  not  participate  in 
fluid  deformation  and  acts  only  to  rigidly  rotate  a  fluid 
element.  The  effect  of  this  term  defines  the  fluid  vor- 
ticity  and  can  be  shown  to  be  antisymmetric. 

rdv.  dV-'. 

^  =  V2  d  -  (Ai  .26a) 
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(A1.26b) 


The  antisymmetry  of  the  vorticitv  tensor  insures 
that  only  three  independent  components  exist.  It  is  some¬ 
times  helpful  to  define  a  vorticity  vector  such  that: 


Z.  =  e  .  . .  uj .  . 
k  13k  1 j 
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1  i^j^k  ;  i,j,k  form  a  right-handed 
system . seguence 

=  -f-1  i^j^k  ;  i,j,k  form  a  left-handed 
^  system.  .  .  -anticyclic  seguence 


[  0  everywhere  else .  acyclic  seguence 


Al.1.6  The  Stream  Function 

Consider  the  two-dimensional  flow  field  defined 
by  Figure  Al.8.  The  region  of  interest  for  this  discus- 
sion  lies  between  two  streamlines.  Since  by  definition, 
the  local  velocity  vectors  are  tangential  to  a  streamline, 
no  fluid  particles  may  cross  that  imaginary  boundary. 


Figure  Al.8  Flow  between  two  streamlines 
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Concentrating  on  the  control  volume  ABC,  one  may  apply 
the  conservation  of  mass  principle  for  an  incompressible 
fluid- 

mass  efflux  +  mass  efflux  _  „ 
across  AB  across  BC 

(pv2dx1)  +  (pv1dx2)  =  0  (A1.27) 

The  stream  function  is  defined  according  to  eq.(A1.28). 

dip  =  -v2dx1  =  v1dx2  (A1.2  8) 

Qualitatively  then,  the  stream  function  can  be  viewed 
as  a  measure  of  the  flow  between  two  given  streamlines 
in  the  velocity  field. 


Often  it  is  advantageous  to  define  the  velocities 
in  terms  of  a  stream  function.  It  may  be  possible  in  cer¬ 
tain  problems  to  substitute  the  equivalent  of  the  velocity 
components  in  terms  of  the  stream  function. 


v 


1 


dip 

3x2 


(A1.29) 


This  maneuver  reduces  the  number  of  dependent  variables 
to  be  solved  for  at  the  expense  of  raising  the  order  of 
the  differential  equation.  Despite  the  drawbacks,  many 
problems  have  been  solved  using  this  technique. 
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A1 . 2  Constitutive  Relationships 

The  general  conservation  laws  that  were  developed 
in  the  previous  section  made  no  mention  of  any  particular 
material  system.  Clearly,  when  one  attempts  to  solve 
real  engineering  problems  in  fluid  mechanics,  the  behav¬ 
ior  of  the  fluid  is  an  important  consideration.  One  usu¬ 
ally  wishes  to  relate  the  response  of  the  fluid  to  the 
applied  stresses.  The  expressions  that  link  the  stresses 
and  the  fluid  motion  are  called  constitutive  relationships. 
It  should  be  pointed  out  that  these  relationships  are  not 
derivable  from  basic  conservation  principles,  but  are 
merely  approximations  at  best.  Many  times  the  expressions 
are  such  that  they  seem  to  fit  the  data  over  a  wide  range 
of  deformation  rates.  Other  approximations  result  from 
experience  or  intuition  that  have  worked  well  over  the 
years . 


Al.2.1  The  Inviscid  Fluid 

The  simplest  constitutive  relationship  that  could 
be  proposed  is  one  in  which  the  stress  tensor  is  isotropic, 
namely 

a .  .  =  -p 6 . .  (A1 . 30a) 

xj  ^  ID 

-p  0  01 

=  0  -p  0 

.0  0  -pJ 


(Al. 30b) 
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From  eq. (A1.30)  it  is  apparent  that  the  stress  tensor  is 
independent  of  the  deformation  rate  of  the  fluid.  Also, 
notice  that  the  viscous  effects  do  not  enter  into  the 
rheological  characterization  of  the  fluid,  as  implied  in 
the  definition  of  inviscid. 

If  the  fluid  is  gaseous  and  compressible  then  the 
pressure  may  be  related  to  the  temperature  and  the  density 
via  the  ideal  gas  law 

P/P  =  RT  (A1.31) 

where  R  is  the  gas  constant. 

A  much  broader  class  of  fluids  may  be  investigated 
if  one  assumes  that  the  fluid  approximates  the  incompres¬ 
sibility  criterion.  As  usual,  one  begins  with  the  basic 
conservation  of  mass  principle.  In  general  vector  notation, 
recall  eq. (A1.7)  given  by 

V* (pq)  =  0  (Al. 32) 

where  V  is  defined  to  be  the  del  operator  and  q  is  the 
velocity  vector.  Incompressibility  states  that 

p  =  constant 

and  so  eq.CA1.32)  becomes 


V  •  q  =  0 


(Al.  33) 
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This  last  equation  states  that  the  flow  field  is  solenoi- 
dal  and  it  is  seen  that  the  divergence  of  the  velocity 
vector  is  zero  everywhere. 

Now  if  the  vorticity  (defined  by  eq.(A1.26))  van¬ 
ishes  identically  throughout  the  flow,  it  will  be  said 
that  the  flow  field  is  irrotational .  In  fact,  if  the  vor¬ 
ticity  is  initially  zero  for  our  inviscid  fluid,  even  then 
the  velocity  field  is  irrotational.  To  see  this,  suppose 
that  initially  (either  in  time  or  space)  the  flow  possesses 
no  rigid  body  rotation.  The  only  way  that  a  fluid  element 
could  assume  angular  motion  is  through  the  action  of  sur¬ 
face  tractions  in  the  form  of  shear  stresses  along  its 
edges.  Clearly  the  assumption  of  inviscid  flow  prohibits 
the  development  of  any  such  shear  stresses.  This  fact  was 
borne  out  by  the  constitutive  relationship  of  eq.(A1.30). 
Mathematically,  for  a  flow  field  to  be  irrotational,  the 
curl  of  the  velocity  field  must  vanish. 

Vxq  =  0  .  (A1.34) 

In  addition,  eq.(A1.34)  assures  us  that  a  velocity  poten¬ 
tial  exists  such  that 

q  =  V$  .  (A1.35) 

Inserting  eq.(A1.35)  into  eq.(A1.33)  yields  the  well-known 
Laplace  equation  which  governs  potential  flow  theory. 
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Notice  the  simplifications  that  result  from  the  above 
substitution.  The  governing  differential  equation  con¬ 
tains  only  the  unknown  scalar  function  3>  instead  of  three 
unknown  velocity  components , 

=  0  .  (A1.36) 

Comparing  with  eq.(A1.32),  one  sees  that  a  trade-off 
exists  in  that  the  order  of  the  differential  equation 
has  been  increased.  Fortunately,  many  solution  techniques 
are  available  for  solving  the  Laplace  equation. 

The  addition  of  appropriate  boundary  conditions 
establishes  the  well-posed  nature  of  such  problems  solved 
under  these  assumptions.  Potential  theory  is  applicable 
to  flows  in  which  the  fluid  viscosity  is  negligible.  A 
convenient  non-dimensional  parameter  defined  to  be  the 
Reynold's  number 

Re  =  inertial  forces 
viscous  forces 

determines  the  range  of  applicability  of  potential  theory. 
Large  values  of  the  Reynold's  number  are  indicative  of 
flow  conditions  dominated  by  inertial  effects  and  thus  the 
approximation  of  an  inviscid  fluid  is  quite  good. 


10.3 


A1 .2.2  The  Newtonian  Fluid 

Unlike  the  inviscid  fluid  of  section  A1 . 2 . 1 ,  the 
Newtonian  fluid  behaves  as  a  viscous  fluid.  The  stress 
tensor  is  linearly  related  to  the  rate  of  deformation  by 
a  single  constant  known  as  the  viscosity  of  the  fluid. 

a  .  =  -P6.  .  +  2yd.  .  (A1.38) 

When  one  limits  the  response  of  the  stress  tensor  to  a 
linear  dependence  on  the  rate  of  deformation  tensor,  the 
claim  is  made  that  the  parameter  y  is  not  a  function  of 
d^ j .  The  viscosity  may  still  depend  upon  the  temperature, 
however.  Quite  often  the  temperature  effect  upon  the  vis¬ 
cosity  may  be  an  important  consideration  and  not  at  all 
negligible. 

Returning  to  eq.(A1.38),  one  may  note  that  a  scale 
factor  of  2  is  placed  in  front  of  the  second  term.  This 
is  purely  for  convenience  in  that  it  exactly  cancels  the 
l/z  which  is  found  in  the  definition  of  d_^  (eg.  (A1 . 26)  )  . 
When  the  Newtonian  constitutive  relationship  is  substi¬ 
tuted  into  the  equations  of  motion  the  resulting  expres¬ 
sions  are  known  as  the  Navier-Stokes  equations.  They  are 
quite  important  and  have  received  widespread  attention 
over  the  years.  Often  they  are  used  as  a  first  approxi¬ 
mation  to  difficult  flow  problems  as  a  means  of  obtaining 
useful  insight. 
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A1 . 2 . 3  The  Generalized  Newtonian  Fluid 


Here,  the  viscosity  is  seen  to  depend  upon  the 
deformation  rate  tensor  as  indicated  by  eq.(Al.39), 

Gij  =  _p6ij  +  ndij  (A1.39) 

where  r)  =  variable  viscosity  function. 


We  anticipate  expressing  the  viscosity  function  n  in 
terms  of  the  invariants  of  the  deformation  rate  tensor. 


It  is  a  logical  choice  since  one  expects  the  constitutive 
relation  of  eq.(A1.39)  to  be  unchanged  by  a  coordinate 


transformation . 
shown  below. 


The  three  invariants  of  the  tensor  d  are 

% 


Hd  =  V2  (  (tr  d)  2  -  tr(d)2) 

IHfl  =  det(d) 

'v* 


(A1.40) 

(A1.41) 

(A1.42) 


where  "tr"  is  the  abbreviation  for  the  trace  of  a  matrix. 
Simply,  it  is  seen  to  be  the  sum  of  the  components  along 
the  main  diagonal.  Since  it  has  been  decided  that  the 
viscosity  should  be  some  function  of  the  above  invariants, 
one  proceeds  to  determine  this  dependence. 


Looking  at  the  first  invariant,  one  easily  recog¬ 
nizes  this  to  be  the  conservation  of  mass  condition.  For 
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Cartesian  coordinates  and  assuming  the  fluid  to  be  in¬ 
compressible,  this  was  shown  to  be  identically  zero  in 
eq.(A1.8).  As  a  result,  it  is  permissible  to  ignore  any 
dependence  from  Id  on  the  viscosity.  It  is  also  custom- 
ary  to  ignore  the  IIId  dependence  upon  the  viscosity. 

'V 

In  doing  so  one  is  neglecting  the  extensional  flow  effects 
upon  the  rheological  approximation.  Essentially  then,  the 
model  will  be  limited  to  shearing  type  flows.  This  is  an 
acceptable  simplification  in  lieu  of  the  fact  that  many 
important  problems  may  still  be  solved  if  extensional  ef¬ 
fects  are  ignored. 

Thus  far,  it  has  been  shown  that  the  viscosity 

function  for  the  generalized  Newtonian  fluid  depends  on 

the  second  invariant  of  the  deformation  rate  tensor  II,  . 

d 

n  =  f(Hd)  (A1.43) 

'V- 

There  have  been  many  models  which  have  specified  the 
exact  nature  of  eq.(A1.43)  in  terms  of  the  flow  conditions. 
Most  of  the  models  are  empirical  in  nature  and  have  been 
derived  as  a  result  of  curve  fitting  schemes  with  adjust¬ 
able  parameters.  The  best  expressions  have  a  minimum  of 
parameters  that  must  be  determined  experimentally  and  also 
exhibit  a  wide  range  of  applicability.  Perhaps  the  best 
known  two-parameter  model  is  that  of  the  power  law  fluid. 
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The  viscosity  function  for  the  power  law  model  is 


given  below  for  laminar  shearing  flow. 


n  =  M| T |n  1  (A1.44) 


The  constant  M  corresponds  most  naturally  to  the  viscosity 
for  the  Newtonian  fluid,  while  T  represents  the  shear  rate 


of  the  fluid.  The  shear  stress  is  related  to  an  arbitrary 


power  of  the  deformation  rate,  hence  the  nomenclature: 


power  law  model.  Comparing  this  expression  to  eq.(A1.38) 
it  is  evident  that  for  n  =  1  the  equations  are  identical 


and  the  Newtonian  constitutive  relation  is  recovered. 


Figure  A1.9  shows  the  variation  of  shear  stress  with  de¬ 


formation  rate  for  both  Newtonian  and  power  law  fluids. 


Figure  A1 . 9  Variation  of  shear  stress  with 

deformation  rate  for  a  power  law 
fluid 
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Values  of  the  power  law  index  (n)  less  than  unity 
correspond  to  pseudoplastic  or  shear  thinning  behavior. 
Index  values  greater  than  one  represent  dilitant  or  shear 
thickening  fluids. 

An  example  is  given  below  which  illustrates  the 
derivation  of  the  power  law  model  for  steady-shearing 
flow.  Consider  the  two-dimensional  Couette  type  flow 
shown  in  Figure  A1.10.  The  kinematics  of  the  flow  are 
quite  straightforward  and  are  listed  here  for  convenience. 

v,  =  f(^2)  (A1.45a) 

v?  =  0  .  (A1.45b) 


Figure  A1.10  Two-dimensional  Couette  flow 
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The  only  deviatoric  stress  resulting  from  this  type  of 
flow  happens  to  be  To  justify  this  claim,  recall 

that  from  eq.(A1.39) 

P. .  =  nd. . 

1J  ij 

and  that  the  only  nonzero  component  of  the  deformation 
rate  tensor  for  the  above  kinematics  was  d12- 

f  ,  1 

o  y2r  0 

<*ij  =  Vzi  0  0  (A1.46) 

0  0  0 

l 

From  eq.(A1.41),  it  is  easy  to  calculate  the  second 
invariant  of  the  deformation  rate  tensor 

JId  =  dlld22  "  d12 

— h-Hr* 

and  by  solving  for  the  shear  rate  (T)  ,  one  finds 

T  =  2/-IEd  (A1.47) 

'V 

Substituting  eq.(A1.47)  into  eq.(A1.44),  the  viscosity 
function  is  now  found  in  terms  of  the  second  scalar 
invariant. 

n-1 

_  ,nn-l,  _  .  2 

M2  ( tEE  d ) 

v 


0 


(A1.4  8) 
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Again  note  that  if  n=  1,  n  —  M  =  constant  then  the  viscos¬ 
ity  function  reduces  to  the  value  of  the  Newtonian  vis¬ 
cosity. 


Thus,  one  may  summarize  the  findings  of  the  above 
discussion.  As  written,  the  power  law  model  is  generally 
applicable  to  the  shearing-type  flows,  hence  stress  terms 
such  as  P.^  ,  P22  ,  etc.  do  not  arise.  Secondly,  for  vis¬ 
cous  flows  that  satisfy  the  kinematics  of  eq.(A1.45)  the 
rheological  constitutive  relationship  becomes 

p.  .  =  Mir |n 
ID  1  1 

Al.2.4  The  Oldroyd- Maxwell  Fluid 

So  far  the  constitutive  equations  that  have  been 
discussed  have  been  independent  of  time.  The  stress  ten¬ 
sor  has  been  shown  to  depend  upon  the  deformation  rate 
tensor  for  both  Newtonian  and  power  law  fluids.  Now,  the 
stresses  are  allowed  to  depend  upon  the  time  domain. 
Physically,  the  stresses  will  be  time-dependent  if  the 
fluid  possesses  some  elasticity.  In  this  case  one  would 
expect  normal  deviatoric  stresses  to  be  present  and  thus 
give  rise  to  anomalous  effects. 

In  essence,  a  relationship  is  proposed  to  describe 
time- dependent  rheological  behavior  which  suggests  the  idea 


of  viscoelasticity.  In  fact,  the  Maxwell  model  is  prob¬ 
ably  the  best  known  viscoelastic  model  and  also  the  easi¬ 
est  to  understand.  The  constitutive  relationship  for  the 


0 1 dr oyd- Maxwell  fluid  is  given  by  eq.(A1.49 


P.  .  +  e: 


2yd.  . 

1 3 


(A1.49) 


With  the  exception  of  the  second  term,  eq.(A1.49)  is 
identical  to  that  of  the  Newtonian  fluid.  In  eq.(A1.49) 

0  is  defined  to  be  a  relaxation  time  constant.  It  may  be 
further  defined  as  the  viscosity  divided  by  the  shear 
modulus  (y/G)  which  has  the  units  of  time.  The  derivative 
Wt  is  defined  as  the  Oldroyd  convective  derivative  (OCD) . 
Assuming  Cartesian  coordinates,  the  OCD  may  be  expanded  to 


Ap  .  .  3P  .  . 

Ill  =  _ hi 

At  3t 


9P  •  .  3v .  3v. 

+  v  -  p .  — 1  -  p  — - 

m  3x^  lm  3x  mi  3x 

m  m  J  m 


.  (A1.50) 


As  before,  it  is  convenient  to  describe  the  stresses  in  a 
moving  coordinate  system.  The  variations  in  the  stress 
tensor  with  respect  to  time  are  then  given  by  the  OCD. 

The  first  term  of  the  OCD  (eq.(A1.50))  records  the  vari¬ 
ation  in  stress  with  respect  to  the  local  coordinate  sys¬ 
tem.  The  remaining  terms  take  into  account  the  time  vari¬ 
ations  as  a  result  of  the  fluid  velocity  and  nonlinearities, 
respectively.  Together,  they  model  the  total  derivative 
of  stress  with  respect  to  time. 


Ill 


To  best  illustrate  the  implications  of  the  Maxwell 
model  (eq . (Al . 49 ) ) ,  an  example  is  worked  out  below.  Con¬ 
sider  the  same  two-dimensional  steady,  shearing,  flow  of 
Figure  A1.10.  Recall  that  the  kinematics  were 


V1  =  rx2 
v2  =  0 


and  the  deformation  rate  tensor  was  shown  to  be 


o  y2r 
y2r  o 


0  0 


Both  the  kinematics  and  the  strain  rate  tensor  are  iden¬ 
tical  for  this  problem.  The  difference  arises  when  one 
applies  a  constitutive  relationship.  Applying  eq.(A1.49] 
with  the  known  components  of  d  ,  it  is  not  hard  to  show 
that  the  only  nonzero  stress  terms  are  contained  in  the 
equations 


Pn  +  -  0 

&12 

P12  +  e_^T~  " 


(Al . 51a) 


(Al ,51b) 


Working  out  the  Oldroyd  derivatives ,  one  sees  that 


P—  =  -2P  - — 

ml  3x 

m 


-2 


(Al . 52a) 
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If  one  compares  the  structure  of  the  stress  ten¬ 
sor  with  that  of  the  deformation  rate  tensor,  an 

interesting  observation  can  be  made.  The  normal  stress 
P^1  is  nonzero  even  though  the  d^  component  of  the  d 
tensor  is  zero.  What  is  happening  is  that  the  Maxwell 
fluid  develops  inherent  normal  stresses  when  subjected  to 
pure  shearing  flow.  This  anomalous  effect  is  common  among 
many  polymeric  solutions.  In  many  instances,  these  normal 
stresses  may  be  far  greater  than  the  shearing  stresses 
needed  to  deform  the  material;  hence,  they  may  not  be 
neglected. 


113 


A1 . 3  Selected  Analytical  Examples 


In  an  effort  to  reinforce  some  of  the  concepts 
presented  in  this  discourse  and  also  illustrate  the 
method  of  problem  solving,  two  examples  are  presented 
in  this  section.  These  problems  serve  to  introduce  the 
reader  to  the  basic  idea  of  simple  shearing  flow  in  a 
channel.  Contrasting  behavior  in  the  fluid  velocity 
fields  result  from  employing  two  different  constitutive 
equations . 

As  one  shall  see,  solving  problems  in  fluid 
mechanics  is  basically  intuitive.  However,  there  is  a 
logical  sequence  of  events  that  should  be  discussed 
briefly.  First,  one  must  make  an  educated  guess  as  to 
the  flow  kinematics.  This  initial  guess  can  be  only  as 
good  as  one's  insight  or  experience  may  dictate.  Clearly, 
the  kinematics  that  are  chosen  must  satisfy  the  boundary 
conditions  of  the  physical  problem.  The  next  step  usu¬ 
ally  involves  choosing  a  suitable  constitutive  equation. 

If  this  is  not  explicitly  given  in  the  problem  definition, 
one  must  rely  on  past  experience.  With  the  assumed  fluid 
behavior  in  hand,  the  above  information  is  substituted 
into  the  equations  of  motion.  This  step  acts  to  solve 
for  the  missing  information  in  terms  of  unknown  stresses 
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and/or  the  exact  expressions  for  the  velocity  components. 
If  the  initial  guess  as  to  the  kinematics  was  incorrect, 
that  fact  will  frequently  be  borne  out  during  the  solu¬ 
tion  to  the  equations  of  motion.  Since  the  details  of 
the  flow  problem  are  a  direct  consequence  of  the  assumed 
kinematics,  incorrect  guesses  as  to  the  kinematics  result 
in  invalid  and  often  contradictory  details.  If  need  be, 
added  information  from  the  continuity  equation  can  be 
used  as  a  check  to  verify  the  final  solution.  One  shall 
see  how  to  apply  this  methodology  in  the  problems  that 
follow. 

Al.3.1  Planar  Poiseuille  Flow 

Consider  the  flow  of  a  very  viscous  fluid  between 
two  stationary  walls  as  shown  in  Figure  Al.ll. 


vrf  (|-‘52  ) 


Figure  Al.ll  Normalized  velocity  profile  for 
two-dimensional  Poiseuille  flow 
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As  instructed,  a  set  of  reasonable  kinematics  is  proposed. 

V1  =  (A1 . 54a) 

v2  =  0  (A1.54b) 

v3  =  0  (A1.54c) 

In  the  problem  formulation  it  is  assumed  that  the  dimen¬ 
sion  of  is  sufficiently  large  so  that  no  variation  in 
the  flow  takes  place  in  that  direction.  In  fact,  one 
further  specifies  that  in  addition  to  being  a  constant, 
is  zero.  For  clarity,  one  many  draw  on  the  analogy 
from  the  plane  strain  approximation  in  linear  two- 
dimensional  elasticity.  The  velocity  component  in  the 
x2  direction  is  also  seen  to  be  zero  since  one  expects 

no  velocity  fluctuations  transverse  to  the  flow  direction. 
Finally,  v^  is  considered  to  be  a  function  of  x2  since  it 
is  known  from  the  "no  slip"  boundary  condition  at  the 
solid  wall  that  v^  =  0  at  this  interface  and  yet  v^  must 
be  finite  within  the  channel  for  flow  to  take  place.  The 
velocity  component  v^^  is  not  a  function  of  x^  since  one 
is  assuming  shearing  flow.  Hence,  it  may  be  justified 
that  the  kinematics  of  eqs.(A1.54)  are  reasonable. 

The  equations  of  motion  for  steady  inertialess 
flow  in  Cartesian  coordinates  are  written  below. 
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(A1 . 55a) 


(A1.55b) 


(A1.55c) 


Usually,  when  the  transient  and  inertial  terms  (the  left- 
hand  side  of  the  equations  of  motion)  are  ignored,  the 
flow  is  such  that  viscous  effects  dominate  the  inertial 
effects.  Recall  from  section  A1 . 2 . 1  the  definition  of 
the  Reynold's  number.  For  slow,  viscous  flows  the  Re 
number  approaches  zero  and  the  equations  of  motion  reduce 
to  those  of  eqs.(A1.55). 


From  the  kinematics,  eqs.(A1.54),  one  may  at  once 
write  down  the  deformation  rate  tensor. 


3  v .  3  v  . 

dij  =  1/2  [sJJ  +  3^ 


0  0 


(A1.56) 


0  0  OJ 


Upon  employing  the  Newtonian  constitutive  relationship 

o .  .  =  -p<5 .  .  +  2yd .  . 
i3  *13  H  ID 


(A1.57) 
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it  may  be  concluded  that  the  only  nonzero  deviatoric 
stress  component  is  is  due  to  the  linear 

relationship  between  the  P  and  d  tensors.  Hence,  eqs . 

%  "u 

(A1.55)  become 


0 


-dp  ,  9 

dx^  9x2 


3v. 


0 


( 2  y  Vz 


dp 

dx2 


(A1.58a) 


(A1.58b) 


respectively.  Eq.(A1.55c)  provides  no  information  in  that 
it  is  identically  zero;  thus  it  will  be  ignored  in  the 
remaining  discussion. 


From  eq.(A1.58b)  it  is  apparent  that 

dvl 

dx„  “  f(x2) 


(A1.59; 


which  reinforces  the  original  choice  of  kinematics.  To 
solve  for  v^  eq.(A1.58a)  is  integrated  twice  with  respect 


to  x. 


x?  C, 


r  =1JP_3  +  Iix  +12 

1  p  dx1  2  \i  2  v 


(A1.60) 


Upon  the  application  of  the  two  boundary  conditions 


v^  =  0 
dv. 


dx. 


x_  =  ±a 


x2 


(A1 . 61a) 


0 


0 


(A1.61b) 
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it  turns  out  that  c1 = 0  and  c2 =  (-a  /2)dp/dx1  where 
"a"  is  defined  as  half  of  the  channel  width.  Therefore, 
eq.(A1.60)  becomes 


v  =  -L  -dP_  (x2 
1  2y  dxx  X2 


(A1.62) 


Since  (dp/dx.)  is  negative  (a  pressure  drop  is  exper¬ 
ienced  for  viscous  flow),  one  may  multiply  eq.(A1.62) 
by  (-1)  and  deal  with  |  dp/dx.jJ  which  may  be  written 
simply  as  dp/dx^  with  no  confusion  intended.  Upon  re¬ 


arranging  eq . (A1 .62) 
a  ^  do 

vi  =  2S  11 


(A1.63) 


If  one  considers  a  flow  with  a  centerline  velocity 
of  unity,  then  setting  x2 = 0  in  eq.(A1.63)  it  is  evident 
that  the  collection  of  constants 

2  j 

a  dp 

2y  dx^ 

must  equal  1.  Now  one  may  use  the  remaining  continuity 
equation  to  determine  the  average  value  of  the  velocity. 


v.dx_  =  v 
1  2  avg 


(A1.64) 


so  that  upon  carrying  out  the  integration  one  finds 
Vava  =  2/3  vmax  *  At  last'  the  desired  expression  for 
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the  normalized  velocity  profile  is  obtained  for  Newtonian 


channel  flow. 

2 

vx  =  3/2  (1  ~  -§  )  (A1.6  5) 

& 


This  result  illustrates  the  well-known  parabolic  velocity 


profile  characterisitc  of  planar  Poiseuille  flow.  This 


relationship  is  referenced  several  times  throughout  the 


preceding  chapters. 
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r>  =  Ml 


Recall  that  M  is  interpreted  to  be  the  consistency  of  the 
fluid  and  n  is  defined  to  be  the  power  law  index.  Sub¬ 
stituting  for  P^2  in  sq .  (A1.66)  yields  a  second  order 
ordinary  differential  equation  in  terms  of  v. 


jsl  .  utLfUr' 

dx^  dx2  [dx2J 


(A1.68) 


Integrating  eq. (A1.68)  twice  and  making  use  of  the  bound¬ 
ary  conditions  of  eqs.(A1.6i)  the  resulting  velocity 


field  is  shown  to  be 


n+1  n+1 


v  =  _2 _ 

1  n+1  1/n  lax 


1J  ^ 


(A1.69) 


Resorting  to  similar  arguments  as  in  the  previous  section, 
let  it  be  assumed  that  the  maximum  centerline  velocity  is 
equal  to  unity.  The  average  velocity  may  be  calculated 
from  the  continuity  equation  to  be 


v  =  (1 - - —  )  v 

avg  2n+l  ;  max 


(A1.70) 


As  before,  one  may  normalize  eq.(A1.69)  with  respect  to 
the  average  velocity  to  obtain  the  general  velocity 
profile  for  a  power  law  fluid. 


l+2n  2 

V  -|  =  -T- -  1 - 

1  1+n  a  ; 


(A1.71) 
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In  general,  these  velocity  profiles  will  not 
be  parabolic  except  in  the  special  case  where  n  =  1  and 
eq.(A1.71)  collapses  to  that  derived  for  the  Newtonian 
fluid  (eq. (A1 . 65) )  .  Figure  Al.12  shows  the  variation 
of  the  velocity  profile  with  power  law  index.  These 
are  the  characteristic  profiles  for  pseudoplastic  or 
shear  thinning  fluids.  We  will  be  concerned  primarily 
with  the  response  of  such  fluids  undergoing  shearing 
deformation. 


0  x,  1.0  2.0 


Figure  Al.12  Velocity  profile  for  the 
power  law  fluid 
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FINITE  ELEMENTS  IN  FLUID  MECHANICS 


A2.1  Introduction 


One  important  aspect  of  the  determination  of 
fiber  orientation  within  a  flow  field  relies  on  the 
ability  to  completely  solve  the  fluid  mechanics  problem. 

By  completely,  it  is  meant  that  all  the  components  of 
the  velocity  field  (v. )  and  the  independent  components 
of  stress  (a^j)  are  known,  as  well  as  approximations 
to  the  stream  function  and  vorticity.  In  Appendix  Al, 
it  was  apparent  that  only  a  limited  number  of  problems 
could  be  solved  by  obtaining  analytical  solutions  for 
the  velocity  field.  For  a  simple  two-dimensional  flow 
field  the  equations  of  motion  were  integrated  directly. 

As  the  kinematic  expressions  increase  in  complexity, 
such  a  simplified  approach  fails.  Even  for  Stoke' s  flow, 
where  the  equations  of  motion  are  relatively  simple,  it 
becomes  formidable  to  obtain  a  closed  form  solution  if 
the  domain  geometry  becomes  too  complex.  As  one  antici¬ 
pates  the  nonuniform  mold  shapes  that  may  arise,  it  becomes 
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necessary  to  seek  an  approximation  to  the  actual  solution. 
The  finite  element  method  shall  be  employed  for  this  pur¬ 
pose.  Irregular  boundaries  will  be  shown  to  present  no 
difficulty  with  this  type  of  formulation.  In  fact,  through 
proper  discretization  some  rather  exotic  boundary  shapes 
may  be  modeled  quite  closely,  thus  insuring  a  good  repre¬ 
sentation  of  the  solution  domain. 

In  discussing  the  application  of  the  FEM  to 
specific  problems,  one  shall  assume  that  the  reader 
has  a  limited  but  working  knowledge  of  the  subject. 

Many  excellent  texts  are  listed  in  the  references  for 
introductory  or  additional  reading.  What  is  intended 
here  is  to  present  a  logical  sequence  of  steps  that  will 
provide  a  basis  toward  understanding  the  FEM  as  applied 
to  fluid  mechanics  for  this  work.  Many  techniques  such 
as  domain  discretization,  the  assembly  of  elements,  etc., 
are  direct  carry-overs  from  the  solid  mechanics  area. 

Other  similarities  will  be  mentioned  as  they  arise  to 
assist  the  reader  in  his  understanding.  The  main  differ¬ 
ence  between  solid  and  fluid  mechanics  formulations  lies 
in  the  derivation  of  the  element  equations.  Again,  the 
material  contained  in  Appendix  A1  will  be  of  paramount 
importance  in  the  development  of  such  equations. 


124 


In  section  A2 . 2 ,  a  quick  glance  at  the  varia¬ 
tional  methods  is  presented  in  relation  to  the  finite 
element  method.  It  is  useful  to  introduce  variational 
methods  to  present  an  alternate  approach  to  the  general 
formulation  via  the  direct  stiffness  approach.  Most 
treatises  introduce  the  FEM  by  the  direct  stiffness 
approach  because  of  its  simplicity.  Our  goal 
is  to  compare  the  formulations  due  to  the  direct  stiff¬ 
ness  approach  and  variational  methods.  Also  an  illumi¬ 
nating  example  is  included  to  illustrate  the  link  between 
the  governing  differential  equation  and  the  discretized 
matrix  representation  of  a  problem. 

The  six  node  triangular  element  is  developed  in 
section  A2.3.  Included  within  this  scope  are  references 
to  the  standard  element  and  the  interpolation  functions. 
Having  chosen  a  suitable  element  with  prescribed  shape 
functions,  the  actual  solving  of  the  differential  equa¬ 
tion  is  anticipated.  The  method  of  weighted  residuals 
(MWR)  allows  one  to  approximate  the  solution  within  the 
domain  to  an  artibrary  degree  of  accuracy.  In  section 
A2 . 4 ,  a  particular  form  of  the  MWR,  the  Galerkin  method, 
is  introduced  for  this  purpose. 
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The  element  equations  governing  Stoke' s  flow 
are  derived  in  section  A2.5.  Ideas  from  sections  A2 . 3 
and  A2 . 4  are  utilized  to  arrive  at  the  final  set  of 
discretized  equations.  The  solution  of  the  linear 
system  of  equations  is  enhanced  with  a  frontal  elimina¬ 
tion  method,  and  the  procedure  is  outlined  in  section 
A2.6.  A  procedure  that  utilizes  the  Galerkin  method 
and  determines  the  behavior  of  the  derived  quantities 
of  stream  function  and  vorticity  is  outlined  in  section 
A2.7.  These  two  functions  are  needed  in  the  development 
of  the  fiber  orientation  model.  Lastly,  an  example  is 
solved  via  the  FEM  that  draws  on  many  of  the  ideas  pre¬ 
sented  in  Appendix  Al.  A  simple,  two-dimensional, 
shearing  flow  problem  is  solved  numerically  to  show 
the  validity  in  the  FEM. 

A2 . 2  The  Finite  Element  Approach 


A2.2.1  Direct  Stiffness  Approach 

There  are  two  possible  avenues  to  traverse  when 
introducing  the  finite  element  method.  The  first,  and 
perhaps  the  most  widely  adopted,  is  called  the  direct 
stiffness  approach.  Its  derivation  originated  in  the 
area  of  solid  mechanics  primarily  because  initial 
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investigations  into  the  finite  element  approach  were 
made  in  this  area.  However,  direct  analogies  exist  in 
other  disciplines  as  well,  which  should  emphasize  that 
it  is  a  rather  general  methodology.  This  procedure 
shall  be  briefly  outlined  below,  for  the  purpose  of 
comparing  this  method  with  that  of  the  variational 
approach  to  be  discussed  subsequently.  As  one  will 
appreciate,  each  approach  has  its  advantages,  but  both 
paths  of  pursuit  yield  identical  results. 

Consider  the  one-dimensional  domain  D  shown  in 
Figure  A2 . 1 .  For  definiteness  it  may  be  considered  to 
be  a  uniform  rod  subject  to  compressive  forces.  The 
independent  variable  is  defined  as  x  while  the  dependent 
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Figure  A2 . 1  One- dimensional  domain,  D, 
for  the  FEM 
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or  solution  variable  shall  be  denoted  as  y.  Since  it 
was  chosen  to  consider  a  rod  in  compression,  the  vari¬ 
able  y  many  be  viewed  as  in-plane  displacements.  To 
illustrate  the  discretization  procedure  one  may  divide 
the  original  domain  D  into  four  elements  that  connect 
a  total  of  five  nodes.  Each  node  in  this  simple  exam¬ 
ple  is  to  possess  one  degree  of  freedom  (i.e.  a  single 
displacement) .  A  typical  element  is  shown  in  Figure 
A2.2.  Note  that  each  element  possesses  two  degrees  of 
freedom  y^  and  and  that  the  variation  of  the  depen¬ 

dent  variable  is  assumed  to  be  linear  within  the  indiv¬ 
idual  element.  We  shall  expand  upon  the  ideas  presented 
here  later  when  discussing  the  representation  of  the 
unknown  field  variable  in  terms  of  shape  functions. 


y 


Figure  A2.2  A  one-dimensional  element  with 
two  degrees  of  freedom 
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If  one  turns  to  solid  mechanics  to  borrow  a 
constitutive  relation  in  which  to  relate  these  dependent 
variables,  one  finds  from  Hooke's  law  that 

ky  =  f  (A2.1) 

where  k  =  stiffness  coefficient 
y  =  displacement 
f  =  force 

The  analogy  with  forces  and  displacements  is  one  that 
usually  allows  the  reader  to  think  in  terms  of  familiar 
quantities.  Now,  typical  with  the  finite  element  pro¬ 
cess  one  asks  that  eq.(A2.1)  be  satisfied  within  each 
subdivision  or  element.  Thus,  writing  eq.(A2.1)  for 
each  element,  one  arrives  at  a  system  of  simultaneous 
equations  to  be  solved.  In  matrix  or  operator  notation 
these  equations  are  seen  to  be 

nxn  nxl  nxl 

[k]  {y}  =  { f  >  (A2.2) 

The  labor  involved  in  this  procedure  results  from  the 
computation  of  the  individual  "stiffness  coefficients". 

It  is  not  hard  to  show  that  this  approach  will  always 
render  the  stiffness  matrix  symmetric.  An  additional 
important  point  to  realize  from  eq.(A2.2)  is  that  the 
maximum  size  of  the  matrices  depends  upon  the  number 
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of  variables  contained  in  the  overall  system.  In  this 
example,  the  number  of  variables  equals  the  number  of 
nodes . 

A2.2.2  Variational  Methods 

One  begins  this  section  by  considering  a  family 
of  functions  on  the  interval  between  two  given  points, 
say  xq  and  x^  as  shown  in  Figure  A2.3.  Clearly  there 
are  a  vast  number  of  possible  relationships  that  satisfy 
the  above  criterion.  The  number  of  admissible  functions 
may  be  restricted  by  requiring  that  they  all  contain 
continuous  derivatives  of  at  least  second  order.  In 
addition,  it  may  be  assumed  that  they  all  satisfy  the 
essential  boundary  conditions  at  the  two  endpoints. 

Y(x0)  =  yQ  (A2.3a) 

y(xi)  =  yi  (A2.3b) 


Figure  A2 . 3  The  variation  of  y(x)  on 
the  interval  (x^  -  xq) 
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Assume  that  the  family  of  curves  may  be  represented  by 
the  integral 


Ky(x)  ) 


F  ( x,y , y  1  )  dx 


(A2.4) 


where  the  prime  denotes  differentiation  with  respect 
to  x.  The  value  of  the  integral  inherently  depends 
upon  the  particular  choice  of  the  function  F.  The  func¬ 
tion  F  is  seen  to  depend  only  on  x,  y  and  y' .  More  gen¬ 
eral  derivations  include  dependence  upon  higher  order 
derivatives,  but  for  current  purposes  the  noted  depen¬ 
dence  will  suffice.  Also,  the  integral  I  contains  only 
one  independent  variable  whereas  more  general  procedures 
might  include  more.  Ultimately,  I(y(x))  is  seen  to  be 
a  function  of  a  function  and  one  many  define  this  type 
of  relationship  by  the  term  functional. 


It  will  be  of  concern  to  determine  the  stationary 
value  of  the  functional  I.  It  will  be  apparent  that  this 
extremal  condition  will  lead  to  a  differential  equation 
that  must  be  satisfied  if  the  functional  is  to  be  an 
extremum.  This  differential  equation,  often  called  the 
Euler  equation,  will  often  be  identical  to  the  equations 
that  govern  various  physical  phenomena.  Thus,  one  shall 
see  that  for  certain  problems,  it  will  be  possible  to 
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work  with  the  minimizing  of  a  functional  instead  of 
dealing  directly  with  the  differential  equation.  This 
is  not  always  possible  since  differential  equations  that 
govern  some  physical  circumstances  possess  no  variational 
counterpart.  With  the  motivation  now  in  place,  one  may 
therefore  proceed  to  minimize  the  function  of  eq.(A2.4). 

Assume  that  y(x)  is  the  true  solution  in  that 
it  extremizes  I.  Next,  one  allows  for  the  following 
perturbation  on  the  function  y(x)  of  the  form 

y(x)  =  y(x)  +  £n.(x)  (A2.5) 

where  e  is  a  small  number.  The  new  function  n (x) ,  which 
may  be  chosen  arbitrarily,  satisfies  the  boundary  con¬ 
ditions  n(xQ)  =0  and  n  (x-^)  =0.  Schematically,  the  func¬ 
tions  of  eq .  (,A2 . 5 )  are  shown  in  Figure  A2.4.  We  shall 


Figure  A2 . 4  The  function  n (x) 


132 


consider  these  functions  as  fixed  and  the  only  variable 
is  the  scalar  parameter  e.  Substitute  eq.(A2.5)  into 
eq.(A2.4)  and  arrive  at  the  following  function. 

I1 

I  (y  (x)  )  =  F  (x,y+er)  (x)  ,y  ' +en  '  (x)  )  dx  (A2.6) 

x 

o 

Since  eq.(A2.6)  is  a  function  only  of  the  parameter  z, 
one  may  find  the  extremum  by  ordinary  means.  The  deriva¬ 
tive  of  the  function  with  respect  to  the  parameter  e  will 
provide  the  necessary  condition  for  an  extremum  to  occur. 


dl (y  (x)  ) 

ds  e=0 


(A2.7) 


If  one  applies  eq.(A2.7)  to  eq.(A2.6), 


dl (y (x) ) 

d£  £=0 


3F  ]  f_3F 

3y  £=ojr|  3y '  £=0 


D  '  dx 


(A2. 8) 


and  integrating  the  second  term  of  eq.(A2.8)  by  parts 


The  first  term  of  eq.(A2.9)  vanishes  because  of  the 


133 


essential  boundary  conditions  placed  upon  the  function 
r)  (x)  .  The  remaining  terms  must  then  equal  zero  per 
eq . ( A2 . 7 ) . 


'  9F 

,9y 


X 


0 


_d_f_9F 
dx [ 3y ' 


'l 

ri  (x)  dx 


j 


0 


(A2.10) 


The  function  n (x)  is  not  necessarily  zero  everywhere; 
thus  in  order  to  satisfy  the  equality,  the  bracketed  por¬ 
tion  of  eq.(A2.10)  must  be  identically  zero. 


3F  d_f9F  ) 

9y  ~  dx  [  3y 1  J 


(A2.ll) 


Eq.(A2.11)  is  the  Euler  equation  and  is  a  necessary  con¬ 
dition  for  extremizing  the  functional  of  eq.(A2.6).  The 
question  of  a  sufficient  condition  for  minimum  extrema 

may  be  investigated  by  computing  the  second  variation  of 

d2 1 

the  functional  with  respect  to  e,  namely  - — . 

de  *■ 


How  can  the  above  mathematics  be  applied  to  the 
finite  element  method  as  introduced  in  the  previous  sec¬ 
tion?  The  so-called  variational  approach  will  be  illus¬ 
trated  through  the  following  example.  Suppose  one  knows 
that  the  governing  differential  equation  to  a  particular 
problem  admits  the  functional 
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xi 

f  1  7 

I(y(x))  =  j  ( y2  (dy/dx)  -  f  (x)y)  dx  .  (A2.12) 

x 

o 

As  shown  above,  it  is  desired  to  minimize  this  functional 
subject  to  the  appropriate  boundary  conditions.  In  doing 
so,  one  will  satisfy  the  necessary  condition  of  eq.(A2.11) 
and  thus  satisfy  the  differential  equation.  Consider, 
for  definiteness,  the  same  domain  as  outlined  in  Figure 
A2.1.  Since  one  is  about  to  apply  the  FEM,  the  initial 
consideration  is  concerned  with  representing  the  solution 
in  a  discretized  manner.  The  unknown  variable,  y,  may  be 
represented  by  the  series  expansion 

2 

Y  =  I  W.y.  (A2.13) 

i=l  1  1 


within  each  element,  where  y^  denotes  nodal  values  of  the 
unknown  field  variable  and  VA  corresponds  to  the  interpola¬ 
tion  functions.  The  interpolation  functions  dictate  the 
variation  of  the  field  variable  within  any  element.  Assum¬ 
ing  a  one-dimensional,  linear  element  as  described  by 
Figure  A2.2,  the  shape  functions  become 


xi  "  x 

W.  =  — -  =  1  -  l 

1  xi  -  xo 

X  -  X 

w,  =  - - 

2  x .  -  x 


(A2 . 14a) 


=  € 


(A2 . 14b) 
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Figure  A2 . 5  Linear  shape  functions 

Characteristic  of  all  shape  functions,  the  i-th  shape 
function  takes  on  the  value  of  unity  at  the  i-th  local 
node  and  equals  zero  at  all  other  element  nodes.  In 
this  fashion  one  is  able  to  approximate  the  solution  in 
terms  of  the  nodal  values  of  the  unknown  function. 

Also,  the  standard  element  is  introduced  as  a 
means  of  calculating  derivatives  with  repect  to  space 
coordinates.  On  a  given  element  the  normalized  coordi¬ 
nate  £  ranges  from  0  to  1  and  is  related  to  the  global 
coordinate  x  by  the  transformation  shown  in  Figure  A2 . 6 . 
The  advantage  of  working  on  a  standard  element  will  be 
shown  as  opposed  to  working  with  global  coordinate  values. 
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Figure  A2 . 6  Transformation  to  local 
coordinates 


Integrations  which  need  to  be  performed  will  be  over  a 
unit  interval  and  this  idea  is  useful  when  implementing 
numerical  integration  schemes. 


Thus,  returning  to  the  functional  to  be  mini¬ 
mized,  eq.(A2.12)  becomes 


I  = 


(  Vz  (dy/dx)  -  f(x)y)dx  ,  y(0)=y(L)=0  (A2.15) 


where  L  ranges  over  the  full  domain.  Next,  it  is  natural 
to  discretize  the  domain  and  sum  up  the  contributions  from 


all  the  elements. 
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Substituting  eqs.(A2.19)  and  eq.(A2.18)  into  the  functional 
expression  of  eq.(A2.17)  yields 
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x  =  (yi+ryi) 

j  2(x.+1-x.) 


-  f(x(C)) [yi(i-?)+yi+1C] 


:xi+1-x.)d5 


(A2.20) 


or  equivalently 


snr^Ti  (Yi  '  2YiYi+i  +  yi+i: 


-  yi(xi+1-xi)  (!-£;)  f(x(C))dC 


yi+l(xi+l“xi)  f (x(?))Cd5 


(A2.21) 


If  one  defines  the  following  quantities: 


k.  .  =  - - - r 

i,r  (xi+1-x±) 


(A2 ,22a) 


k.  ... 
l ,  i+l 


(xi+l~xi; 


(A2.22b) 


^i+l,i+l  ^i,i 


(A2 .22c) 


Fi  =  (xi+1"xi)  (l-5)f(x(C))dC 


(A2.22d) 


Fi+1  =  (xi+rxi}  J5f  (X(C)  )d? 


(A2.22e) 


then  one  may  write  eq.(A2.21)  in  a  more  convenient  form  as 
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I  . 
J 


'MYi,Yi+1} 


k .  .  k .  >  - 

1,1  1 , l+l 


ki+l,i  ki+l,i+l 


i+1 


■  ‘a-a+i1 


i+ij 


(A2.23) 


Introducing  operator  or  matrix  notation  eq.(A2.23)  can 
be  compacted  to  yield 

Ij  =  1//2yjkjyj  ~  YjFj  ^n°  SUm)  (A2.24) 

where  superscript  T  represents  the  transpose  of  a  matrix. 

To  complete  the  exercise,  it  is  necessary  to  sum  the  con¬ 
tributions  from  eq.(A2.24)  and  then  carry  out  the  minimi¬ 
zation.  Mathematically,  this  correponds  to 


dZl . 
_ I  = 

dyj 


N 


,  y  (  y2y?k  -y  •  “  yXF  . )  =  0 
dy^  [  /2Jj  jJj  ■} 


(A2.25) 


If  it  is  true  that  the  order  of  summation  and 
differentiation  may  be  interchanged  then  eq.(A2.25)  pro¬ 
vides  a  set  of  simultaneous  equations  that  must  be  solved. 


N 


I  ^ — ■(  :/2yTk  .y  .  -  y?F  .  )  =  0 

l  dy^  /2Jj  yj  31 


(A2.26) 


Eq.CA2.261  reduces  to 

nxn  nxl  nxl 
[kj  {y}  =  {F} 


(A2.27) 
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after  carrying  out  the  minimization. 

Clearly,  then,  eq. (A2.27)  is  exactly  the  same  as 
eq.(A2.2)  and  one  concludes  that  the  results  of  the  vari¬ 
ational  approach  are  identical  to  the  direct  stiffness 
approach,  as  is  always  the  case  for  a  linear  system. 

A2.3  Choice  of  Elements 


The  finite  element  method  centers  around  the  abil¬ 
ity  to  find  an  approximate  solution  within  the  solution 
domain.  Typically,  one  divides  the  domain  into  smaller 
divisions  of  a  simpler  shape  and  works  exclusively  to 
approximate  the  solution  within  the  element.  Let  us  brief¬ 
ly  investigate  the  notion  of  interpolation  or  shape  func¬ 
tions  to  determine  how  one  approximates  the  behavior  of 
the  unknown  within  an  element. 

Consider  the  two-dimensional  triangular  element 
shown  in  Figure  A2.7.  Its  simplicity  is  rivaled  only  by 
the  four  node  rectangular  element,  but  it  has  the  distinct 
advantage  of  being  able  to  better  approximate  irregular 
boundaries.  Therefore,  only  the  triangular  element  will 
be  considered  in  this  discussion. 


X 


Figure  A2.7  Linear  triangular  element 
in  global  coordinates 


The  interpolation  functions  are  functions  of  the  coor¬ 
dinates  and  possess  the  inherent  property  that  W.  (x.)=  5.  . 

i  3  11 

This  is  to  say  that  in  eq.  (A2.28)  the  function  <j>  must 
equal  the  field  variable  at  that  node  and  may  not  have 
any  contribution  from  other  nodal  values . 


One  seeks  to  approximate  the  unknown  variation  of 
the  field  variable  <p  throughout  the  triangular  domain  as 

4>  =  Wi<Di  (A2.28) 

where  are  defined  as  the  interpolation  functions  and 
<jK  are  the  nodal  values  of  <j> .  This  will  be  possible  if 
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in  the  case  of  the  triangular  elements  we  use  a  two- 
dimensional  polynomial  to  expand  <j> .  One  retains  only 


<}>  =  a,  +  a_x  +  a,y 


(A2.29) 


the  linear  terms  of  the  expansion  which  admit  three 
unknown  constants:  ,  a2  ,  and  a^.  These  constants 

may  be  exactly  solved  for  since  eq.(A2.29)  applies  to 
each  of  the  three  nodal  values  of  the  element,  namely. 


1  xi  yj  K 


*2  =  1  x2  y2  a2 


1 


(A2.30) 


where 


1  xi  yi 
1  x2  y2 
1  x3  y3 


[CJ  (x.,y.)  =  coordinate 
i  1 i  values 


Upon  inversion  of  eq.(A2.30)  we  find 


{a}  =  [C]  1{4>  > 


(A2 . 31) 


and  one  may  rewrite  eq.(A2.29)  as 


<p  =  [1 , x,y ]  {a}  =  [l,x,y]  [C ] —  x } 


(A2.32) 


If  one  compares  eq. (A2.32)  with  eq. (A2.28)  it  is  apparent 
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that  the  interpolation  functions  are 

W±  =  [l,x,y] [C]"1  (A2.33) 

The  method  outlined  above  for  obtaining  the  interpolation 
functions  is  straightforward  in  principle  but  may  fail 
to  yield  meaningful  information  if  [C]  ^  fails  to  exist. 

We  endeavor  to  reveal  a  method  that  allows  one  to  write 
down  interpolation  functions  by  inspection. 

As  one  has  seen  from  the  previous  section,  it  is 
beneficial  to  define  a  standard  element  whose  normalized 
coordinate  values  range  from  0  to  1 .  The  standard  tri¬ 
angular  element  is  displayed  in  Figure  A2 . 8  with  the 
corresponding  transformation  equations  given  in  deter- 
minental  form  by  eqs.(A2.34), 


n 


Figure  A2 . 8  Standard  triangular  element 
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x-x^ 

X3"X1 

Y-Yf 

Y3’Yl 

X2“X1 

X3_X1 

Y2-Yi 

Y3“Yi 

X2'xl 

x-x^ 

Y2“Yi 

y~yl 

2A  | 


(A2.34a) 


(A2.34b) 


The  denominator  of  eq. (A2.34b)  corresponds  to  twice  the 
area  of  the  triangular  element. 

In  terms  of  these  local  coordinates,  one  seeks 
the  interpolation  functions  that  describe  the  variation 
of  <p  within  the  element.  Consider  and  the  figure  below. 


n 


Figure  A2 . 9  Determination  of  shape  functions 
by  inspection 
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By  definition 

Wx  =  1  ( C=0 ,  n=0)  (A2 ,35a) 

and  = 0  at  points  (0,1)  and  (1,0).  The  second  criterion 
may  be  satisfied  if  one  asks  that  =  0  along  the  line 

f(?,n)  =?+n-l=0  .  (A2.35b) 

The  interpolation  function  that  satisfies  eqs.(A2.35a)  and 
(A2.35b)  becomes 

Wx  =  1  -  £  -  n  •  (A2.36a> 

The  second  interpolation  function  is  simply 

W2  =  K  (A2.36b) 

and  it  can  easily  be  verified  that  W2  =  0  at  nodes  1  and  3 
and  is  unity  at  node  2.  Also,  it  may  be  shown  that 

W3  =  n  -  (A2.36c) 

Often  these  shape  functions  are  called  area  coordinates 
which  satisfy  the  relationships 

Ai 

W]_  =  —  =  l-  ?-  n  (A2.37a) 

A2 

W2  =  —  =  K  (A2 . 37b) 

A3 

=  —  =  r|  (A2.37c) 

where  A  defines  the  area  of  a  triangle.  Clearly  as  a 
consequence  of  eqs . (A2.37)  the  shape  functions  must 
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satisfy  the  relationship 
Wx  +  W2  -f-  w3  =  1 

n 


Figure  A2.10  Area  coordinates 

The  three-node  triangle  has  been  shown  to  model 
the  linear  variation  of  the  field  variable  within  the 
element.  With  this  type  of  approximation  the  field 
variable  is  continuous  at  the  element  interfaces .  One 
often  refers  to  this  type  of  behavior  as  Co  continuity. 

It  is  possible  to  obtain  quadratic  interpolation 
for  the  triangular  element  if  mid-side  nodes  are  added  to 
the  existing  nodes  of  Figure  A2.ll.  The  transformation 
equations  to  the  standard  element  remain  identical  with 
those  of  eqs.(A2.34);  however,  one  would  expect  the  inter¬ 
polation  functions  to  be  different. 


Figure  A2.ll  Six  node  (quadratic)  triangular 
element  in  local  coordinates 

To  see  how  one  may  construct  interpolation  functions 
by  inspection,  consider  the  development  of  for  the 
six  node  triangle. 


n 


Figure  A2.12  Determining  iK  by  inspection 
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The  shape  function  iK  must  vanish  for  nodes  that  lie 
along  line  ab  and  cd.  To  insure  this  fact  one  forms 
a  product  of  these  criteria,  namely 


1>1  *  (1  -  K  ~  n)  (  Vz  ~ 
v  (Wx)  (Wx  -  Vz) 


(A2. 38) 


Patching  up  the  proportionality  with  a  constant  renders 
the  shape  function  1/^  =  1  at  £  =  0 ,  n  =  0  resulting  in  the 
final  expression 


^  =  2W1(W1  -  Vz) 


(A2. 39a) 


The  derivation  of  the  remaining  shape  functions  follow 
similar  arguments  and  are  listed  below. 


iP2  =  2W2(W2  -  l/z) 

ip3  =  2W3(W3  -  Vz) 

iP4  =  4WxW2 

*5  =  4W2W3 
h  =  4W3W1 


(A2 . 39b) 
( A2 . 39  c ) 
( A2 . 39d) 
(A2 . 39e) 
(A2 . 39  f ) 


It  is  logical  to  wonder  what  these  shape  functions  look 
like  that  have  been  derived.  The  shape  functions  iK  and 
ip4  for  the  six  node  triangle  are  shown  in  Figure  A2.13. 

One  may  be  disillusioned  at  this  point  into 
thinking  that  placing  more  nodes  along  element  boundaries 


will  necessarily  increase  the  degree  of  approximation. 


Figure  A2.13  The  shape  functions  and 

for  the  six  node  triangle 

To  explain  why  this  need  not  be  true  we  construct  Pascal 
triangle,  shown  in  Figure  A2.14. 


Figure  A2.14  Pascal's  triangle  illustrating 
the  degree  of  interpolation 
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The  triangular  array  of  variables  lists  the  terms  of  a 
two-dimensional  polynomial  expansion.  The  most  accurate 
results  to  the  FEM  are  obtained  when  the  interpolation 
functions  contain  the  highest,  complete  order  polynomials. 
Inclusion  of  a  few  terms  of  higher  order  in  the  interpo¬ 
lation  functions  tends  not  to  increase  the  accuracy  of 
the  solution  and  often  imparts  noise  into  the  system  of 
equations.  For  example,  suppose  for  certain  problems 
it  is  advantageous  to  work  on  a  four  node  quadrilateral 
element.  The  interpolation  functions  would  contain  the 
terms  l,x,y,xy.  This  type  of  representation  for  the 
field  variable  only  admits  linear  variation  along  element 
sides.  The  inclusion  of  the  xy  term  adds  little  in  the 
form  of  increased  accuracy  to  the  solution  approximation. 

For  the  six-node  triangular  element  Pascal's  tri¬ 
angle  reveals  complete  quadratic  representation  of  the 
field  variable  throughout  the  element.  This  feature  of 
quadratic  interpolation  functions  as  well  as  ease  of 
boundary  modeling  make  the  six  node  triangle  a  desirable 
element  to  use.  This  element  shall  be  employed  when 
solving  the  fluid  mechanics  problem  by  finite  element 


methods . 
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A2 . 4  Method  of  Weighted  Residuals 


We  have  made  reference  to  alternate  methods  of 
solving  a  problem  that  is  governed  by  a  differential 
equation.  In  section  A2.2.2,  it  was  assumed  that  a 
variational  principle  existed  for  the  phenomenon  under 
consideration  and  that  one  could  deal  with  the  minimiza¬ 
tion  of  a  functional.  It  was  also  pointed  out  that  many 
problems  are  not  amenable  to  a  variational  approach  and 
one  is  forced  to  work  directly  with  the  differential 
equation.  In  this  section,  the  method  of  weighted 
residuals  is  illustrated  as  it  applies  to  the  finite 
element  method. 

Consider  the  following  equation 

Lu (x)  =  f  (A2.40) 

where  L  is  an  ordinary  differential  operator.  The  boundary 
conditions  are  such  that 

u(x)  =  u  (A2.41) 

for  all  x  on  the  boundary.  The  possible  functions  u(x) 
under  consideration  must  be  square- integrable  in  the  sense 
that 

o 

2 

|u(x)  I  dx  <  °°  .  (A2 .42) 

a 
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Recall  our  discrete  approximation  to  the  continuous 
function 


u* (x)  =  >  w.u. 

i=i  1  1 


CA2.43) 


In  general,  eq. (A2.43)  will  not  satisfy  the  gov¬ 
erning  equation  exactly.  One  may  anticipate  some  error 
or  residual  to  arise  if  the  approximation  is  substituted 
into  the  governing  differential  equation.  This  may  be 
shown  by  writing  eq.(A2.40)  with  u(x)  being  replaced 


by  u* (x) , 


m*  (x)  -  f  =  R 


/  u . LW .  -  f 


(A2. 44) 


Now,  the  method  hinges  on  whether  one  can  minimize  the 
residual  that  is  found  in  eq. (A2.44) .  To  accomplish 
this  a  set  of  weighting  functions  v^  are  introduced  and 
it  is  required  that  the  weighted  integral  of  the  resi¬ 
dual  vanish,  namely 


<v^,R>  =  v^Rdx  =  0 


(A2.45) 


If  one  imposes  this  condition  on  eq.(A2.44)  the  result 


becomes 


<  [  u. (LW. , v . ) >  -  <f ,v.>  = 

2  _  T  -L  J-  I  | 


0 


(A2. 46) 
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To  this  point  the  weighting  functions  Vj  have 
been  completely  arbitrary  and,  in  fact,  to  specify  them 
will  define  a  particular  form  of  the  method  of  weighted 
residuals.  There  are  several  possible  choices  each 
having  their  own  merit,  but  the  most  useful  seems  to 
be  the  Galerkin  method.  Here,  the  weighting  functions 
are  designated  to  be  the  same  shape  or  interpolation 
functions  that  are  used  to  describe  the  unknown  variation 
of  the  field  variable  within  the  element.  Therefore 
eq.(A2.46)  becomes 

N 

<  l  u. (LW. ,W.)>  -  <f,W.>  =  0  .  (A2.47) 

i=l  1  1  J  J 

To  understand  why  the  Galerkin  method  works,  one 
may  investigate  the  underlying  assumptions  of  eq.(A2.45). 
The  weighting  functions  which  were  chosen  as  the  inter- 
plation  functions  form  a  complete  set  of  functions.  That 
is  to  say  that  any  given  function  may  be  expanded  as  a 
linear  combination  of  the  shape  functions  which  is  borne 
out  by  eq.(A2.43).  Also  the  shape  functions  may  be  seen 
to  be  linearly  independent  if  one  recalls  how  they  were 
derived.  They  assume  a  value  of  unity  at  the  i-th  node 
and  equal  zero  at  all  other  nodes.  Thus,  any  function 
is  capable  of  being  represented  by  eq.(A2.43)  provided 
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sufficient  terms  are  utilized.  Finally  then,  one  may 
claim  that  the  residual  of  eq.(A2.45)  is  forced  to  zero 
since  it  is  a  function  that  must  be  orthogonal  to  every 
member  of  the  complete  set  of  shape  functions . 

For  a  linear  operator  L,  the  Galerkin  method 
will  always  lead  to  the  same  results  as  those  obtained 
via  the  variational  methods.  Also,  it  is  of  interest 
that  the  Galerkin  method  satisfies  the  governing  differ¬ 
ential  equation  in  the  mean.  This  behavior  contrasts 
that  of  other  techniques  which  approximate  the  pointwise 
convergence  of  the  solution. 

A2 . 5  Derivation  of  Stoke' s  Problem 


Thus  far  we  have  been  concerned  primarily  with 
some  preliminary  ideas  that  are  a  part  of  all  finite 
element  problems.  It  seems  appropriate,  now,  that  one 
becomes  specific  and  deals  with  a  particular  problem  in 
an  attempt  to  unify  the  ideas  presented  so  far  and  at 
the  same  time  draw  on  some  important  concepts  discussed 
in  Appendix  A1 .  The  governing  element  equations  for 
two-dimensional  Stoke' s  flow  will  be  derived  in  a  Carte¬ 
sian  (x,y)  reference  frame. 
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One  may  begin  be  recalling  the  governing  equa¬ 
tions  of  motion  for  fluid  mechanics  systems.  For  Stoke' s 
approximation  these  equations  simplify  to 


-3p 

3x 


+ 


3P 


xx 

3x 


+ 


3P 


_xy 

3y 


+  f 


x 


0 


(A2. 48a) 


-3p 

3y 


+ 


3P 


_xy 

3x 


+ 


f 


y 


o 


(A2.48b) 


and  the  continuity  equation  may  be  written  as 


3u  ,  3v 
3x  3y 


0 


(A2.49) 


Three  primitive  variables  are  anticipated  to 
enter  into  the  formulation:  two  velocity  components  and 
the  pressure.  The  stress  components  will  not  be  con¬ 
sidered  as  primitive  since  they  are  related  directly  to 
the  velocity  components  via  the  Newtonian  constitutive 
relations . 


P 

xx 


P 

yy 

p 

xy 


(A2 . 50a) 
(A2 . 50b) 
(A2.50c) 


It  is  natural  to  propose  the  following  expressions  to 
describe  the  variation  of  the  primitive  variables  within 


the  element. 
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6 

u  =  l  u.’l.  (A2 . 51a) 

3=1  3  3 

6 

v  =  I  v.i p.  (A2 . 51b ) 

j  =  l  J  3 

3 

P  =  l  (A2 . 51c) 

j=l  D  J 

Apparent  immediately  is  the  fact  that  one  may  specify 
quadratic  interpolation  for  the  velocity  components  and 
linear  interpolation  for  the  pressure.  This  should  be 
obvious  from  the  order  of  derivatives  that  appear  in 
eqs . (A2 .48).  The  pressure  is  usually  specified  to  be 
one  order  less  than  the  velocity  components  so  that  over¬ 
constraints  do  not  arise.  A  typical  element  showing  the 
unknowns  is  shown  in  Figure  A2.15;  note  that  there  are 
a  total  of  15  degrees  of  freedom  per  element. 


Figure  A2.15  Typical  element  for  the 
u,v,p  formulation 
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If  one  applies  the  Galerkin  method  to  eqs.(A2.48) , 
the  following  expressions  result. 


h 


A'-P  +  pxx»  +  Why  +  fxJ 


3  (~P  +  P  )  +  +  f 


3y 


yy 


9x  xy  y 


dS  =  0 


dS  =  0 


(A2.52a) 

(A2 . 52b) 


*i 


,|H  +  15  Us  =  0 
(.  3x  9yJ 


(A2.52c) 


These  equations  may  be  conveniently  written  as 


3  [il).  (-p  +  P  )  ]  +  ~(i(>.P  ) 

Lri  *  xx'  3y  yi  xy 


3^ 


ax 


3xi(-P  +  Pxx} 


+  Tb.f 


dy  xy  'i  xj 


dS  =  0 


(A2.53a) 


9^. 


TT-[^-(-P  +  P  )]  +  TT-(^.P  ) - ^(-P  +  P  ) 

9y  ^  yy  9x  i  xy  9y  r  yy 


3<J>< 


-  P  +i{>.f  dS=0 

^  r  T  n  t  r 


9x  xy  i  y 


(A2.53b) 


+  i 


■  |S  +  U  dS  =  0 

{ 3x  dyj 


(A2. 53c) 


If  one  concentrates  on  the  first  two  terms  of  eqs.(A2.53a) 
and  (A2.53b),  it  is  apparent  that  they  represent  a  divergence 
of  the  vector  g  defined  as 
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g  =  ^i((-p  +  Pxx)i  +  (Pxv)5 


(A2 . 54) 


Recall  from  the  two-dimensional  form  of  the  divergence 


theorem  that 


V*g  dS  =  n*g  dC 


(A2 . 55 ' 


Here  n  is  defined  as  the  unit  outward  normal  vector. 


If  eq . (A2 . 55)  is  applied  to  eqs.(A2.53)  the  result 


is  given  by 


^  9'^  C 

<-53T ''p+Pxx>  +  <TF'Pxy>  =  <fx'V+!'¥xdc  (A2-56i 


dip  i  dip.  . 

<^3T'Pxy>  +<^T'-P+Pyy>  =  +  J*itdC  (A2 . 56b ) 


. ,  3u  3v^  A 
+  37  =  0 


(A2 . 56c) 


The  brackets  indicate  a  scalar  product  of  the  quantities 

separated  by  commas.  The  notation  t  and  t  represents 

x  y 

the  x  and  y  components  of  the  surface  traction  vector,  an 

idea  borrowed  from  continuum  mechanics.  It  is  well  known 

that  any  combination  of  surface  stresses  may  be  resolved 

into  a  single  traction  vector  whose  components  shall  be 

represented  by  t  and  t  . 

x  y 


The  result  of  the  application  of  the  divergence 
theorem  to  the  momentum  equations  has  been  to  remove  the 
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derivatives  from  the  primitive  variables  and  onto  the 
shape  functions.  Derivatives  of  the  shape  functions  are 
trivial  and  may  be  easily  accessed  during  numerical 
calculations . 

final  step  in  the  derivation  of  the  governing 
element  equations  is  to  implement  the  constitutive  relations 
of  eqs . (A2.50) ,  as  well  as  the  discretized  representation 
of  the  primitive  variables.  These  substitutions  render 
the  governing  equations  of  the  following  form. 


3  if;.  3  if;.  3i|>  •  3  if; .  dip. 

ir-2>JTKL>  +  <-§3r<P-35L>luj  + 


3i|; 

f 

'  (<  'V,Pj 

=  <f  ,lp  .  > 
X  rl 

4- 

U.t  dC 

1  1  1  x 

(A2.57a) 

Si jj  .  dlfj  . 

[<“5F'lJ-55I>Iuj  + 

A 

*|h? 

V 

c 

3  if; .  dip. 

dip 

r 

-  JPj 

=  <f  , i p . > 
y  ri 

+ 

J 

ip.t  dC 

1  y 

(A2.57b) 

dip . 

•[<h’  3x>Iuj  - 

dip. 

[<*i'  3y  >lvj 

c 

0 

(A2. 57c) 

For  convenience,  the  bracketed  quantities  are  defined  as 

3i|k  3if;  .  3 if; .  3 if; . 

Aij  =  +  (A2 . 5! 

H*  Hi  3i|>,  dip  . 

Bij  =  + 


(A2 .58a) 


(A2.58b) 
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Since  the  Stoke' s  equations  are  linear,  the  coefficient 
matrix  of  eq.(A2.59)  is  symmetric.  In  fact  the  same 
equation  may  be  derived  from  the  appropriate  variational 
principle.  The  range  of  indices  of  eq.(A2.59)  may  be 
properly  viewed  in  Figure  A2.16  where  the  size  of  the 
submatrices  are  drawn  to  scale  and  labeled. 


Figure  A2.16  System  of  equations  governing 
Stoke 's  flow 
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A2 . 6  Solution  by  Frontal  Elimination 

There  are  several  possibilities  that  exist  to 
solve  the  linear  system  of  equations 

[k]  {y}  =  (f>  .  (A2.60) 

Conceptually,  the  easiest  way  is  to  multiply  eq. (A2.60)  by 
[k]-1  and  obtain 

{y}  =  [kj -1{ f } 

Quite  often  the  [k]  matrix  is  so  large  that  attempts  to 
invert  it  are  cumbersome  if  not  impossible.  An  alternate 
method  is  that  of  Gauss  elimination  in  which  the  [k]  matrix 
is  reduced  to  triangular  form  and  then  solved  by  back  sub¬ 
stitution.  To  assist  in  the  solution  of  the  eq.(A2.60)  by 
Gauss  elimination,  one  may  make  use  of  the  following  help¬ 
ful  ideas.  The  finite  element  method  will  inherently  pro¬ 
duce  many  zero  components  within  the  coefficient  matrix. 

The  positions  of  these  zero  elements  will  be  off  the  main 
diagonal  as  dictated  by  the  nodal  numbering  scheme.  One  may  take 
advantage  of  this  fact  by  retaining  only  a  banded  portion 
of  the  overall  matrix  which  drastically  reduces  the  stor¬ 
age  requirement  for  the  system  of  equations.  Also  if  the 
matrix  is  symmetric,  it  is  necessary  to  store  only  half 
of  the  members.  These  ideas  are  shown  in  Figure  A2.17. 
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Figure  A2.17  Storing  a  banded,  symmetric 

matrix  (x  indicates  a  nonzero 
member) 

The  bandwidth  is  inherently  dependent  upon  the  numbering 
of  the  degrees  of  freedom,  or  more  commonly,  the  global 
nodes.  If  nodes  are  numbered  erratically,  then  large 
bandwidths  may  arise  resulting  in  little  reduction  in 
the  size  of  the  coefficient  matrix  in  core.  Even  when 
we  utilize  the  banded  matrix  it  may  still  prove  too  large 
for  the  existing  computer.  It  is  thus  desirable  for 
smaller  computing  facilities  to  further  reduce  the  size 
of  the  Ik]  matrix.  One  may  introduce  a  method  of  solution 
that  depends  not  upon  the  node  numbering,  but  upon  the 
numbering  of  the  elements.  Invariably,  this  method  will 
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reduce  the  number  of  equations  in  core  by  working  on  a 
very  small  active  matrix. 

The  frontal  elimination  solution  method  was  first 
introduced  by  Irons  [23]  in  1969  as  an  alternative  way 
of  solving  large  systems  of  equations,  especially  those 
of  finite  elements.  The  technique  will  be  briefly  outlined 
here  since  this  was  the  method  used  in  solving  the  fluid 
mechanics  problem. 

Consider  the  finite  element  representation  of  a 
very  long  domain  as  in  Figure  A2.18.  We  number  the  ele¬ 
ments  in  such  a  manner  so  that  their  sequence  indicates 
the  direction  of  the  front.  The  nodes  are  numbered  so  as 
to  provide  for  conveninece  say  in  plotting  the  mesh  or 
some  other  criterion. 


Figure  A2.18  Element  and  nodal  numbering 

scheme  for  frontal  elimination 


Suppose,  for  ease  of  illustration,  that  there  corres¬ 


ponds  only  one  degree  of  freedom  with  each  node,  then  the 
local  "stiffness"  matrices  of  elements  1  and  2  are  as 
follows : 


1 

2 

54 

1 

k(1) 

kl,l 

k(1) 

*1,2 

k(1) 

1,54 

2 

k(1) 

2,2 

k(1) 

*2,54 

54 

symmetric 

k(1) 

54,54 

1 

54 

53 

1 

54 

rk(2) 

Kl,l 

k(2) 

*1,54 

k(2) 

K54 , 54 

k(2)  "1 

Kl,  53 

k(2) 

54,53 

53 

. 

symmetric 

k(2) 

K53 , 53_ 

Element  1 


Element  2 


If  one  assembles  the  contributions  from  elements  1  and  2 
into  the  global  stiffness  matrix  it  becomes  apparent  that 


1  kf^+k.1 


54  symmetric 


2 

53 

54 

(1) 

1,2 

k(2) 

1,53 

k(1)  +k( 
*1,54  1 

(1) 

2,2 

k(1) 

2,54 

k{2) 

53,53 

k(2) 
k54 , 53 

k(1)  +k 
*54,54+* 

(A2.61) 
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Next,  a  working  or  active  matrix  (K  )  is 

AC  I  J[  VE 

defined  to  contain  the  assembled  contributions  to  k. .. 

ID 


kactive 


1 

2 

53 

54 

1 

Kll 

K12 

Kl  3 

K14 

2 

K22 

K23 

K24 

53 

K33 

K34 

54 

symmetric 

K44 

(A2. 62) 


(The  components  KIJ  (I,J=l-4)  are  defined  by 
association  with  those  of  eq.(A2.61).) 


The  active  matrix  will  contain  only  the  equations  asso¬ 
ciated  with  the  variables  on  the  front.  The  variables 
ahead  of  the  front  will  not  yet  be  assembled  while 
those  behind  the  front  are  eliminated.  To  see  how  this 
elimination  occurs,  refer  to  Figure  A2.18.  Upon  assem¬ 
blage  of  the  first  two  elements,  we  note  that  no  other 
elements  contain  node  1.  The  consequence  of  this  is 
reflected  in  the  fact  that  the  first  row  of  the  global 
stiffness  matrix  is  complete.  That  is  to  say,  there  are 
no  more  K-^N  terms  within  the  system  of  equations.  Thus, 
one  may  eliminate  the  first  row  of  eq.(A2.62)  in  the 
usual  fashion,  first  dividing  by  the  pivot  and  subtracting 
this  equation  from  the  remaining  ones  in  the  active  matrix. 
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kactive 


52 


52 

CN 

L 

53 

54 

2 

K22  - 

K12 

Kll 

K23  K13 

K23  Kll 

K24  -  K14 
Kll 

53 

TC33  -  K13 
K33  K11 

K34  - 
J  Kll 

54 

symmetric 

K44-K14 

Kll_ 

(A2.63) 


The  eliminated  equation  may  be  written  on  disk  to  be  re¬ 
called  later  during  the  back  substitution  process.  There¬ 
fore  the  vacant  first  row  of  the  active  matrix  is  ready 
to  receive  an  additional  equation  as  we  assemble  element 

3.  Eq . (A2 .63)  becomes 
52 


52 

2 

53 

54 

k52 , 52 

k52 , 53 

k52, 54 

2 

K22 

K12 

Kll 

VOTJ  K13 

K23  Kll 

K24-K14 

Kll 

53 

kactive 

k53,53 

.  K33  -  K13 

+  K33  Kll 

k53,54 
+  K34  - 

K14 

Kll 

54 

k54 , 54 

symmetric 

+  K44  - 

K14 

Kll 

(A2.64) 

The  subscripted  stiffness  coefficients  refer  to  the  con¬ 
tributions  from  the  local  stiffness  matrix  of  element  3. 
Upon  assemblage  of  the  first  3  elements,  it  is  noted  that 
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in  addition  to  the  variable  at  node  1,  the  variable  at 
node  53  may  also  be  eliminated.  The  procedures  of  as¬ 
semblage  and  elimination  continue  simultaneously  until 
the  entire  assembly  of  elements  is  complete.  The  pro¬ 
found  advantage  of  this  type  of  scheme  is  seen  from  the 
relatively  small  size  of  the  active  stiffness  matrix 
needed  to  be  solved  after  final  assembly.  For  the  prob¬ 
lem  at  hand,  and  according  to  the  particular  numbering 
scheme,  the  largest  that  the  active  matrix  will  be  is 
5x5,  irregardless  of  the  length  of  the  structure.  It 
is  seen  then  that  the  number  of  variables  residing  on 
the  front  determines  the  final  size  of  the  stiffness 
matrix  that  needs  to  be  stored  in  core. 

If  there  exists  a  drawback  to  the  frontal  elim¬ 
ination  scheme,  it  would  be  in  the  extra  work  needed  to 
devise  an  efficient  algorithm  for  keeping  track  of  the 
variable  placement  in  the  active  matrix.  With  the  plac¬ 
ing  of  equations  in  and  out  of  the  active  matrix  it  is 
essential  that  each  variable  is  sent  to  the  proper  position 
during  evaluation.  With  such  a  code,  the  overall  method 
proves  quite  satisfactory  for  implementation  on  small  or 
large  computer  systems. 
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A2 . 7  Calculating  the  Derived  Quantities 

A2 . 7 . 1  Stream  Function 

In  the  previous  section,  the  primitive  variables 
u,  v  and  p  were  solved  for  by  using  the  frontal  elimina¬ 
tion  method.  Now  it  is  desirable  to  know  the  behavior 
of  certain  derived  quantities.  Derived  quantities  shall 
be  defined  as  those  variables  which  are  not  primitive  and 
which  require  additional  numerical  methods  to  discern 
their  behavior.  In  particular,  the  stream  function  and 
vorticity  will  be  referred  to  as  derived  quantities  since 
they  may  be  calculated  in  terms  of  the  primitive  variables. 
We  shall  not  consider  the  stresses  as  requiring  further 
attention  since  their  evaluation  is  seen  directly  through 
constitutive  equations.  What  little  assistence  is  needed 
in  calculating  the  stresses  will  follow  directly  from  the 
ensuing  discussion.  Recall  from  section  Al.1.5  that  the 
velocity  components  were  related  to  the  stream  function 
through  the  equations 

=  ff  and  v  =  ~ff  (A2.65) 

Since  one  has  only  a  discretized  representation  of  the 
velocity  field,  a  method  is  sought  that  will  yield  the 
stream  function  if  the  nodal  velocity  components  are 
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known.  One  such  procedure  utilizes  the  Galerkin  method 
and  is  outlined  below. 


The  starting  point  of  the  derivation  is  to  note 
that  the  Laplacian  of  the  stream  function  may  be  written 
as 


V 


-  9v 
9x 


+ 


9u 

9y 


(A2 .66) 


In  the  usual  manner,  the  approximations  to  the  velocities 
and  stream  function  are  defined  to  be 


6 

u  =  t  u.<j>.  (A2.67a) 

i=l 

6 

v  =  l  v.  4> .  (A2.67b) 

i=l  1  1 

6 

*  =  I  (A2.67C) 

i=l  1  1 


where,  as  before,  are  the  quadratic  interpolation  func¬ 
tions  given  by  eqs.(A2.39).  Now,  substituting  eqs.(A2.67) 
and  also  applying  the  Galerkin  criterion  to  eq.(A2.66),  the 
resulting  expression  becomes 


(A2. 68) 


The  bracketed  notation  indicates  the  appropriate  scalar 
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product  between  the  quantities  separated  by  commas  (i.e. 

6  6 

,  X  E  (p  .  (p  . )  .  In  view  of  the  fact  that  one  is 
1  J  1=13=1  1  J 

dealing  with  discretized  quantities  this  definition  of 
the  scalar  product  follows  directly  from  the  more  famil¬ 
iar  one  of  continuous  quantities  (i.e.  (<p,<p)  E  ff<p<p dxdy)  . 

Removing  the  summation  signs  and  resorting  to 
the  implied  summation  convention,  the  left-hand  side  of 
eq. (A2 . 68)  yields 

'Pj<<Pi,V-V  4>j>  (A2.69a) 

which  may  be  conveniently  rewritten  as 

tp  .  <  V  •  ( (p  .  ,  7d)  .  )  -  7p.,74>.> 

or 

-i|/ j  < V (f> ±  ,  V<j>  j  >  +  <7  •  ( <}>i  ,  7^  j  <p  j  )  >  (A2 . 69b ) 

Employing  the  divergence  theorem  to  eq.(A2.69b),  the 
resulting  equation  becomes 


or 


or 


~'Pj<y<P±r  V(f>  )  *n> 


dip 

j  <7$^ ,  V4>  j  >  +  <  <1)^  > 


f  d,fji 

-q<vt'vq> + 

c 


(A2.69c) 


where  the  underlined  quantities  are  taken  to  be  over  the 
domain  boundary. 
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Upon  substitution  of  eq.(A2.69c)  for  the  left- 
hand  side  of  eq.(A2.68),  we  arrive  at  the  desired  expres¬ 
sion  for  calculating  the  stream  function. 


<  v<#>i ,  v<f>  j  >!f^  =  <<j> 


d<p. 


d(t>. 


i'T i 


>u. 

J 


dip . 


i  3n  Yj 


<P,ac 


(A2.70) 


Immediately,  it  is  evident  that  the  above  expression  is 
of  the  form 


A .  .  ip  . 

i]  ] 


C  . 

l 


where 


(A2.71) 


=  <VcpifV(pj> 
-  A  ^  j 


3d)  . 

<d>  .  ,  -t— 3. 
3y 


>  u . 
D 


f  dlpj 

C 


and  may  be  solved  via  the  methods  of  section  A2.6.  From 
the  above  discussion  it  is  clear  that  the  stream  function 
may  only  be  calculated  subsequent  to  the  solution  of  the 
velocity  field  since  a  knowledge  of  the  nodal  velocity 
components  is  needed  in  the  right-hand  side  of  eq. (A2.71) . 

As  a  final  comment  on  the  general  solution  of  eq.(A2.71), 
the  contour  integration  found  in  eq.(A2.70)  is  identified 
to  be  the  natural  boundary  condition  on  the  stream  func¬ 
tion  and  is  zero  if  the  streamlines  intersect  the  boundaries 
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at  right  angles.  On  portions  of  the  boundary  where  this 
is  not  the  case,  essential  boundary  conditions  (ip  is 
specified)  are  usually  prescribed. 

A2.7.2  Vorticity 

The  vorticity  vector  has  been  previously  defined 
according  to 

<3v.  Sv.'i 

Z,  =  e  .  . ,  co  .  .  =  y2e  .  . ,  i 

k  ljk  ij  i3k(3Xj  3x^J 


In  the  same  manner  as  the  previous  section,  one  seeks  a 
discretized  representation  for  the  vorticity  in  terms  of 
the  primitive  variables.  The  approximating  functions  are 
introduced  in  eqs.(A2.73). 

3 

Z  =  l  Z.’F.  (A2 . 73a) 

1  3  3 

6 

u  =  £  u  .  <p  .  (A2 . 73b ) 

1  3  3 

6 

v  =  £v.<f>.  (A2.73c) 

1  3  3 


Note  here  that  the  representation  of  the  vorticity  vector 
is  linear  in  nature  while  the  usual  quadratic  interpolation 
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functions  are  employed  to  describe  the  velocity  varia¬ 
tions.  This  type  of  formulation  will  yield  nodal  values 
of  the  vorticity  function  at  the  vertices  only.  The 
derivatives  of  the  velocity  components  in  eq.(A2.72) 
suggest  a  lesser  order  of  approximation  to  the  vorticity 
function  to  insure  a  compatible  field.  This  procedure 
is  common  practice  in  numerical  approximation  techniques. 


Substituting  eqs.(A2.73)  into  eq.(A2.72)  and 
applying  the  Galerkin  criterion,  the  following  expres¬ 
sion  results  for  the  vorticity  representation. 


<¥. ,Y  ->Z  . 
i  3  3 


3<fi  . 

=  y2 (u.<¥. , , " 

3  l  3y 


-*->  -  v  .  <  T  . 

j  i  3x 


(A2.74) 


If  one  defines  the  following  quantites 


A.  .  5  <¥.  ,¥  .> 

13  13 

3<j>  .  3 . 

C.  =  y2  (u.<1f.  ~  v.<V. 

l  3  l  3y  j  l  3x 


then  eq.(A2.74)  may  be  written  compactly  as 


(A2. 75a) 
(A2. 75b) 


A.  .  Z  . 
13  3 


and  solved  by  the  methods  described  previously.  With 
the  conclusion  of  this  section,  the  reader  now  possesses 
the  tools  to  interrogate  the  response  of  the  primitive 
variables  as  well  as  the  behavior  of  the  derived  quan¬ 
tities.  In  the  final  section,  we  will  investigate  the 
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solution  to  a  problem  using  the  aforementioned  solution 
methods . 

A2 . 8  Numerical  Solution  of  2-D  Poiseuille  Flow 

Consider  the  shearing  flow  bounded  by  two  station¬ 
ary  walls  as  shown  in  Figure  A2.19a.  The  accompanying 
boundary  conditions  of  the  flow  domain  are  also  illus¬ 
trated  in  the  adjacent  figure.  Note  that  symmetry  has 
been  used  to  reduce  the  area  that  must  be  discretized. 


(b) 

Figure  A2 . 19  (a)  2-D  Poiseuille  flow 

(b)  Boundary  conditions  for 
2-D  Poiseuille  flow 
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One  possible  discretization  is  shown  in  Figure  A2.20a. 

The  domain  consists  of  12  elements  with  quadratic  interpola¬ 
tion  for  both  velocity  components  and  a  linear  interpola¬ 
tion  for  the  pressure.  As  a  result,  the  system  to  be 
solved  contains  82  degrees  of  freedom,  but  using  the  front¬ 
al  elimination  method  the  largest  active  matrix  is  only 
20x20.  The  numerical  results  of  the  velocity  field  agree 
exactly  with  those  imposed  on  the  boundary  as  they  should 
since  the  problem  is  steady.  The  usual  output  form  for 
finite-element  programs  is  of  a  plotted  nature,  which  pro¬ 
vides  a  quick  and  easy  way  to  investigate  the  behavior 
of  the  variable.  Following  this  lead,  contours  of  all 
the  nonzero  quantities  of  interest  for  this  problem  are 
illustrated  in  Figures  A2.20b-f.  As  a  standard,  we  shall 
plot  10  contours  of  the  field  variable,  equally  incre¬ 
mented  between  the  maximum  and  minimum  values.  Having 
viewed  the  behavior  of  the  individual  field  variables  of 
Figures  A2.20b-f,  one  will  find  them  to  be  identical 
with  those  anticipated  from  the  analytical  expressions. 

This  example  provides  the  necessary  confidence  in  the 
solution  of  the  fluid  mechanics  problem  by  the  finite- 


element  method. 


AXIAL  VELOCITY 


SHEAR  STRESS 


Figure  A2 . 20 

Contours  are  plotted  in  10  equal  increments  spanning 
the  range  between  maximum  and  minimum  values . 

(continued  next  page) 


MAX*  0 .0000 


MIN*  -1.0000 


STREAM  FUNCTION 


VQRTICITY 


MIN*  -1.5000 


MAX*  0.0000 


(Figure  A2.20,  continued) 
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C*********************************************************** 

Q*********************************************************** 

c 

c 

C  PROGRAM  "COLINS"  PROVIDES  A  NUMERICAL  SCHEME 

C  THAT  ALLOWS  GRAPICAL  REPRESENTATION  OF  THE 

C  FIBER  ORIENTATION  WITHIN  A  2-DIMENSIONAL 

C  MOLDED  SHORT-FIBER  COMPOSITE. 

C 

C  A  COMPLETE  DESCRIPTION  OF  THE  FLOW  FIELD  MUST 

C  BE  CONTAINED  IN  THE  FILE  "NEWSOL.RES"  PRIOR 

C  TO  IMPLEMENTING  THIS  PROGRAM.  A  FILE 

C  "GRID.DAT"  SUPPLIES  THE  FINITE  ELEMENT  MESH 

C  PATTERN  AND  TOPOLGICAL  INFORMATION. 

C 

C  ALL  NECESSARY  CHANGES  TO  THE  PROGRAM  ARE 

C  DESIGNATED  IN  THE  FOLLOWING  MANNER: 

C 

C  C  $  ***( INSTRUCTION)*** 

C 

Q*********************************************************** 

Q*********************************************************** 


c 

c 

c 

c 

c 

c 


c 

c 


DIMENSION  IC0UNT(20) 

DIMENSION  LTRARR(lO) 

DIMENSION  INDEX(4,3) 

DIMENSION  MENU ( 6 ) 

DIMENSION  XX(3),YY(3),ZZ(3) 

DIMENSION  SXX(3) ,SYY(3) ,SXY(3) ,U(3) ,V(3) ,VORT(3) 


C  $  ***INTRODUCE  THE  NUMBER  OF  NODES (NNODE)  29  TIMES*** 

c************************************************************ 

DIMENSION  XP(99),TP(99),F(99) 

DIMENSION  SXXP(99) ,SYYP(99) ,SXYP(99) ,STTP(99) 
DIMENSION  UP (99) , VP ( 99 ) ,VORTP(99) 

COMMON  /A1 /NNODE , NVAR , NELEM , NBD , NSO , NVEL 
COMMON  /A2/X(99) 

COMMON  /A3/I(99) 

COMMON  /A5/SXXR(20, 99) 

COMMON  /A6/SYYR(20, 99) 

COMMON  /A7 / SXYR (20,99) 

COMMON  /A8/STTR(20,99) 

COMMON  /A9/UR(20,99) 

COMMON  /A10/VR(20,99) 
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COMMON  /A  1 1  ,/VCRTR  (20,99) 

COMMON  /A  1  2  /XPl  (20,99) 

COMMON  /A1 3/YE (20 ,99) 

COMMON  /A 1 4/SXXAE (20,99 ) 

COMMON  /A1 5/SYYAR ( 20 , 99 ) 

COMMON  / A 1 6/SXYAR(20 , 99 ) 

COMMON  /A 1 7 /UAR (20,99 ) 

COMMON  /A1 8/VAR (20, 99) 

COMMON  /A 1 9/V0RTAR(20,99) 

COMMON  /A20/XAR (20,99) 

COMMON  / A2 1 /YAR ( 20 , 99 ) 

COMMON  /A26/PKI ( 1 C , 20 , 60 ) 

C************************************************************ 

C  $  *** INTRODUCE  THE  NUMBER  OF  NODES (NNCDE)  IN  SUBROUTINE  BOUNP 

Q***************4e’*4f’***************************4e-************** 

COMMON  /ICOPY/ICOPY, IFIB 
COMMON  /ITER/ITER ,  OCOUNT 
COMMON  /IHERM/IHERM 

C************************************************************ 

C  $  *** INTRODUCE  THE  HUMBER  OF  ELEMENTS (NELEM)  1  TIME*** 

0************************************************************ 

COMMON  /A4/IN0D(6 , 40) 

Q************************************************************ 

C  $  ***INTR0DUCE  NVAR+NNODE+NNOBE  *** 
0************************************************************ 

COMMON  /D3/Z(1851) 

COMMON  /A23 /BX , BY , MINX , MINY , XMIN , YMIN , SCALE , MENU 
COMON  /A30/HEIGHT 

DATA  ( (INDEX (I , J) ,1=1 ,4),J=1 ,3)/1 ,4,5,4,4,5,3,2,6,6,6,57 

TYPE  1000 

ACCEPT  1010, SC 

TYPE  1020 

TYPE  1030 

ACCEPT  1040,BX,BY 

IF  (BX.EQ.O)  BX=  1. 

IF  (BY.EQ.O)  BY=  1  . 

TYPE  1050 
ACCEPT  1060, ICOPY 
TYPE  1070 
ACCEPT  1060, IAX 

OPEN (UNIT=29 , DEVICE= ' DSK ' ,FILE= ' VARI . DAT ' ) 

0PEN(UNIT=24 ,FILE= ’ GRID.DAT ' ) 

READ (24,1080)  NSO , NNODE , NBD , NVAR , NELEM 

TYPE  1110 

ACCEPT  11 60, IFIB 

TYPE  1120 

ACCEPT  11 60, ITER 

TYPE  1130 

ACCEPT  1 1 60, IHERM 

DO  10  1=  1, NELEM 

10  READ (24, 1 080)  ( INOD( J, I) , J=1 ,6) 
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CONTINUE 

DO  20  1=  1 , NNODE 

20  READ (24, 1090)  X(I),Y(I) 

CONTINUE 

CL0SE(UNIT=24) 

C 

C 

C  CALCULATING  THE  SCALED  PLOTTER  COORDINATES 

C 

C 

DO  30  1=  1, NNODE 

XP(I)=  BX*X ( I ) 

30  YP( I )=  BY*Y ( I ) 

CONTINUE 
XMIN=  XP(1  ) 

XMAX=  XP(1  ) 

YMIN=  YP(  1  ) 

YMAX=  YP( 1 ) 

DO  40  1=  2 , NSO 
XI=  XP(I) 

YI=  YP(I) 

IF  (XI.LT.XMIN)  XMIN=  XI 
IF  (XI.GT.XMAX)  XMAX=  XI 
IF  (YI.LT.YMIN)  YMIN=  YI 
IF  (YI.GT.YMAX)  YMAX=  YI 
40  CONTINUE 

S1=  6.5/(XMAX-XMIN) 

S2=  8 . / ( YMAX-YMIN ) 

IF  (SI . GT.S2)  50,60 
50  SCALE=  S2 

GO  TO  70 
60  SCALE=  SI 

70  ¥IDTH=  SCALE* (XMAX-XMIN) 

HEIGHT =  SCALE*( YMAX-YMIN) 

MINX=  (8.5-WIDTH)/2. 

MINY=  ( 1 0. -HEIGHT )/2.+1 . 

DO  80  1=  1, NNODE 

XP(I)=  MINX+ (XP ( I ) -XMIN )*SCALE 
80  YP(I)=  MINY+(YP(I)-YMIN)*SCALE+1 . 

CONTINUE 

0PEN(UNIT=1 7 , DEVICE* ’ DSK ' ,FILE= ’ NEWSOL. RES ' ) 

C  $  ***FOR  STOKES  PROBLEM  NN=NVAR+5*NN0DE 

C  FOR  MAXLIN  PROBLEM  NN*NVAR+2*NN0DE 

C  FOR  MAXVEL  PROBLEM  NN=NVAR+2*NN0DE 

NN=  NVAR+5*NN0DE 
C 
C 

C  READ  IN  FLOW  DATA  GENERATED  BY  FINITE  ELEMENT  PROGRAM 
C 
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C 

READ (17,11  40 )  (Z(  I ) ,  1  =  1  ,  Nil ) 

CL0SE(UNIT=1 7 ) 

90  TYPE  1 1 50 

C 

C 

C  DETERMINE  THE  FIELD  VARIABLE  FOR  PLOTTING 

C 

C 

ACCEPT  1 1 60, IFUNC 
IF  (IFUNC. EQ. 9)  GO  TO  1 40 
IF  (IFUNC. LE. 4)  100,110 
100  11=  (IFUNC-1 )*NNODE 

GO  TO  150 

110  IF  (IFUNC. EQ. 8)  130,120 

120  11=  (IFUNC-2+IAX)*NN0DE 

GO  TO  150 

130  11=  (5+IAX)*NN0DE+NS0 

GO  TO  150 

140  11=  ( 6 + 1 AX ) *NN 0  DE  +NS  0 

150  111=  NNODE 

IF  (IFUNC. EQ. 7)  111=  NSO 
DO  160  1=  1 ,111 
K=  I 1+1 

160  F( I )=  Z(K) 

CONTINUE 

C 

C 

C  ASSIGN  ALL  NODAL  VARIABLES  FOR  LATER  USE  IN  DETERMINING 
C  INTERPOLATED  VALUES  ALONG  STREAMLINS 

C 
C 

DO  170  1=  1 , NNODE 
SXXP(I)=  Z(I) 

J1 =  NNODE+I 
SYYP(I )=  Z(J1) 

J2=  2*NN0DE+I 
SXYP(l)=  Z(J2) 

J3=  3*NN0DE+I 
UP(I)=  Z(J3) 

J4=  4*NN0DE+I 
VP(I)=  Z(J4) 

J5=  6*NN0DE+NS0+I 
V0RTP(I)=  Z(J5) 

170  CONTINUE 
TYPE  1170 
ACCEPT  11 80, MENU 
TYPE  1190 
ACCEPT  1 200, IDEN 
FMIN=  F( 1 ) 

FMAX=  F  ( 1 ) 
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DO  180  1=  1,NN0DE 
FF=  F( I ) 

IF  (FF.LT.FMIN)  FMIN=  FF 
180  IF  (FF. GT.FMAX)  FMAX=  FF 

CONTINUE 

TYPE  1210, FMIN , FMAX 
TYPE  1220 
ACCEPT  1 270 , IFM 
TYPE  1230 

ACCEPT  1240, FIN, FINCR,NLIN 
TYPE  1250 

ACCEPT  1240, HIN , HINCR , NHLIN 
IF  (ICOPY.EQ. 1 .AND. IFIB.EQ. 1 )  GO  TO  240 
IF  (ICOPY.EQ. 1)  190,200 
190  CALL  PLTSRTC JONES’ ,30) 

GO  TO  210 

200  CALL  PLTSRTC  JONES' ,20) 

210  CALL  FACTOR (SC) 

C 

c 

C  PLOT  MARGIN  AND  BOARDERS 

C 

C 

220  CALL  PL0T(0. ,0. ,3) 

CALL  PL0T(8. 5,0. , 2) 

CALL  PL0T(8.5,11 -5,2) 

CALL  PLCT(0.,11 .5,2) 

CALL  PL0T(0.,0.,2) 

CALL  PL0T(0. ,1 . ,3) 

CALL  PL0T(8. 5 , 1 . ,2) 

CALL  SYMB0L(0. 5, 0.375, 0.225, MENU, 0. ,30) 
C 
C 

C  PLOT  OF  THE  FINITE  ELEMENT  BOUNDARY 

C 

C 

XI =  XP(1 ) 

Y1 =  YP(1 ) 

CALL  PLOT (XI ,Y1 ,3) 

NBD2=  NBD/2 
DO  230  1=  1 , NBD2 
J=  1+1 

IF  (I.EQ.NBD2)  J=  1 
XXX=  XP(J) 

YYY=  YP(J) 

230  CALL  PL0T(XXX, YYY, 2) 

CONTINUE 

C 

C 

C  IDENTIFICATION  OF  CONTOUR  LINES 
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C 

C 

IF  (IFM.EQ.O)  GO  TO  240 

CALL  SYMBOL (-5,1  .25,0.15,  ,0.  ,4) 

CALL  NUMBER (1 .25 , 1 . 25 ,0. 1 5 ,FMIN ,0. , 4) 
CALL  SYME0L(2.75,1  .25,0.15,  ’MAX-’  ,0.  ,4) 
CALL  NUMBER(5*5,1 • 25 ,0. 1 5 ,FMAX, 0. , 4) 

240  CONTINUE 

C 

C 

C  SEARCH  OF  CONTOUR  LINES  BEGINS 

C 

C 

DO  500  IZZ=  1 , NLIN 
JCOUNT-  1 

ZLINE=  FIN+FINCR*( IZZ-1 ) 

DO  290  IELEM=  1 , NELEM 

IF  (IDEN.EQ. 1 )  GO  TO  270 
DO  260  1=  1 ,4 

DO  250  J=  1 ,3 

KK=  INDEX ( I, J) 

L=  INOD(KX, IELEM) 
SXX(J)=  SXXP(L) 

SYY( J)=  SYYF(L) 

SXY( J)=  SXYF(L) 
U(J)=  UP ( L ) 

V(J)=  VP(L) 

VORT( J)=  VORTP(L) 
XX(J)=  XP(L) 

YY ( J ) =  YP(L) 
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ZZ(J)=  F(L) 

CONTINUE 

IZ=  IZZ 

260 

CALL  TRIANG ( ZLINE , XX , YY , ZZ , SXX , SYY , SXY , U , V , 

1  VORT, JCOUNT, IZ) 

CONTINUE 
GO  TO  290 

270  DO  280  1=  1,3 

L=  INOD(I , IELEM) 
SXX(I)=  SXXP(L) 

SYY(I)=  SYYP(L) 

SXY ( I )  =  SXYP(L) 

U(I)=  UP(L) 

V(I)=  VP(L) 

V0RT(I)=  VORTP(L) 
XX(I)=  XP(L) 

YY(I)=  YP(L) 

280  ZZ(I)=  F(L) 

CONTINUE 
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IZ=  IZZ 

CALL  TRIANG ( ZLINE , XX , YY , ZZ , SXX , SXY , SXY , U , V , VORT , 
1  JCOUNT, IZ) 

290  CONTINUE 

ICOUNT(IZ)=  JCOUNT 
300  CONTINUE 


IF  (IRASE.EQ.O)  GO  TO  310 
CALL  ERASE 
CALL  HOME 
310  CONTINUE 

CALL  ERASE 
GO  TO  320 
320  CONTINUE 

C 
C 

C  RESCALE  TO  THE  ORIGINAL  COORDINATES 

C 

C 

DO  330  J=  1 , NLIN 

DO  330  1=  1 ,IC0UNT(J)-1 

XR  ( J , I  )=  ( ( ( XR ( J , I ) -MINX ) /SCALE ) +XMIN ) /EX 
YR ( J , I ) =  ((YR(J,I)-MINY-1 .) /SCALE +YMIN) /BY 
330  CONTINUE 

CONTINUE 

IF  (IC0PY.EQ.1 .0R.IFIB.EQ.1 )  GO  TO  340 
DO  340  J=  1 , NLIN 

WRITE(29> 1 290) 

C 

C 

C  WRITE  DATA  CHECK  IN  FILE  VARI.DAT 

C 

C 

DO  340  1=  1 ,IC0UNT(J)-1 

WRITE(29,1280)  XR( J, I) , YR( J, I) ,SXXR( J, I) , SYYR ( J, I ) , 
1  SXYR (J,I),UR(J,I),VR(J,I), VORTR ( J , I ) 

340  CONTINUE 

CONTINUE 

IF  (ICOPY.EQ. 1 )  GO  TO  350 

OPEN (UNIT=39 , DE VI CE= ' DSK ' , FILE  = ' ORIENT . DAT ' ) 

CL0SE(UNIT=29) 

350  IF  (IFIB.EQ. 1 )  CALL  ARRANG(IC0UNT, NLIN) 

IF  (IFIB.EQ. 1)  CALL  BOUNP(ICOPY,SC) 

IF  (ITER.EQ.O)  GO  TO  360 
DO  360  OCOUNT=  1 , ITER 

IF  (IFIB.EQ. 1)  CALL  INTEG( ICOUNT , NLIN ) 

360  CONTINUE 

CALL  HOME 

IF  (IFIB.EQ. 0)  GO  TO  370 
CALL  BOUNP(ICOPY,SC) 
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CALL  VALUER ( ICOUNT , KLIN) 

IF  (IKESM.EQ.C)  GO  TO  370 

CALL  HERKP  ( NLIN ,  IC CUNT ,  HIS ,  KINCR ,  NHLIN ) 

CALL  HOME 

IF  (IC0PY.EQ.1 )  GO  TO  370 
CLOSE(UNIT=39) 

370  STOP 

1000  FORMAT ( '  INTRODUCE  SCALING  FACTOR') 

1010  FORMAT (F3 . 1 ) 

1020  FORMAT ( '  ENTER  BLOW  UP  FACTOR  IN  X  AND  Y  ') 

1030  FORMAT ( '  **.*  **.*  ’) 

1040  FORMAT (2F5- 2) 

1050  FORMAT ( '  DO  YOU  WANT  HARDC0PY?YES=1 , N0=0 ' ) 

1060  FORMAT(II) 

1070  FORMAT ( '  IS  IT  AN  AXISYMMETRIC  PR0BLEM?Y=1 , N=0 ’ ) 

1080  F0RMAT(6I4) 

1090  FORMAT (2F7. 3) 

1100  FORMAT ( E 1 0 . 3 ) 

1110  FORMAT ( '  DO  YOU  WANT  FIBER  PLOTS?  YES=1,N0=0') 

1120  FORMAT ( '  NUMBER  OF  INITIAL  FIBER  ORIENTATIONS') 

1130  FORMAT ( '  DO  YOU  WANT  A  HERMAN  PL0T?YES=1 , N0=0 ' ) 

1140  F0RMAT(E1 3 • 6 ) 

1150  FORMAT ( 

r  SELECT  F  :SXX=1 ,SYY=2,SXY=3,STT=4,U=5,V=6,P=7,PSI=8,V0RT=9 ' ) 
1160  FORMAT(II) 

1170  FORMAT ( '  INTRODUCE  TITLE: 30  CHARACTERS') 

1180  FORMAT(6(A5) ) 

1190  FORMAT ( '  DO  YOU  WANT  LINEAR (1 )  OR  QUADRATIC (2 )  INTERPOLATION?' 

D 

1200  FORMAT(II) 

1210  FORMAT ( '  FMIN=’ ,E9-3, '  FMAX=',F9.3) 

1220  FORMAT ( '  DO  YOU  WANT  FMIN , FMAX  ON  YOUR  DRAWING?Y=1 , N=0 ' ) 

1230  FORMAT ( '  +0.***E+**+0  .***E+**NN=INI'TIAL  VAL. ,  INCREM .  ,N0  OF  LINES 

D 

1240  FORMAT (2E 10. 3  > 12) 

1 250  FORMAT ( '  +0.***E+**+0.***E+**NN=INITIAL  VAL; , INCR; ,N0  H  LINES' 

1  ) 

1260  FORMAT ( '  CONTOURS  ,E9-3) 

1270  FORMAT (II) 

1280  FORMAT ( 1 X , 8E 1 2 . 4 ) 

1 290  FORMAT(// ///////////////////) 

END 

C 

C 

C 

Q************************************************************ 

c 

c 

C  SUBROUTINES 

C 
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C 

C************************************************************ 
Q******************** **************************************** 

c 

c 

c 

c 

SUBROUTINE  TRIANG ( Z , XX , YY , ZZ , SXX , SYY , SXY , U , V , VORT , JCOUNT , IZ ) 
C 

c 

c 

c 

c********** 

0********** 

C  SCANS  EACH  TRIANGLE  FOR  THE  PROPER  CONTOUR  VALUES  OF 
C  THE  STREAMLINES.  ALSO,  CALCULATES  INTERPOLATED 

C  VALUES  OF  ALL  FIELD  VARIABLES  WHERE  A  STREAMLINE 

C  INTERSECTS  AN  ELEMENT  BOUNDARY. 

q********** 

0********** 

C 

DIMENSION  XX(3),YY(3),ZZ(3) 

DIMENSION  SXX(3) ,SYY(3) ,SXY(3) ,U(3) ,V(3) ,VORT(3) 

COMMON  /A5/SXXR(20,1 ) 

COMMON  /A6/SYYR(20, 1 ) 

COMMON  /A7/SXYR(20,1 ) 

COMMON  /A9/UR(20,1) 

COMMON  /A10/VR(20, 1  ) 

COMMON  /A1 1 /V0RTR(20, 1  ) 

COMMON  /A12/XR(20,1  ) 

COMMON  /A1 3/YR(20 , 1  ) 

COMMON  /ICOPY/ICOPY, IFIB 
IP=  0 

DO  100  M=  1 ,3 
K1=  M 

K2=  M+1-3*(M/3) 

A=  ZZ(K1) 

B=  ZZ(K2) 

C=  (Z-A)*(Z-B) 

IF  (C.GT.O.)  GO  TO  100 

IF  ( Z . EQ . A . AND . Z . NE . B . OR . Z . EQ . B . AND . Z . NE . A )  GO  TO  100 
XI  =  XX (K1  ) 

X2=  XX (K2) 

Y1  =  YY(K1  ) 

Y2=  YY(K2) 

SXX1=  SXX(K1 ) 

SXX2=  SXX(K2) 

SYY1  =  SYY  (K1  ) 

SYY2=  SYY(K2) 

SXY1 =  SXY (K1  ) 

SXY2=  SXY (K2) 
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U1=  U(K1) 

U2=  U(K2) 

VI =  V(K1 ) 

V2=  V  ( K2 ) 

V0RT1 =  V0RT(K1 ) 

V0RT2=  V0RT(K2) 

IF  (A.EQ.E)  GO  TO  50 
X=  XI +(Z-A)*(X2-X1 )/(B-A) 

Y=  Y1 +(Z-A)*(Y2-Y1 )/( B-A) 

SXXR(IZ, JCOUNT)-  SXX1 +(Z-A)*(SXX2-SXX1 )/(B-A) 

SYYR ( IZ , JCOUNT )=  SYY1 +(Z-A)*(SYY2-SYY1 )/(B-A) 
SXYR(IZ, JCOUNT)-  SXY1 +  ( Z-A)*(SXY2-SXY1 )/(B-A) 

UR (IZ, JCOUNT )=  U1 +(Z-A)*(U2-U1 )/(B-A) 

VR(IZ, JCOUNT )=  VI +(Z-A)*(V2-V1 )/( B-A) 

VORTR(IZ , JCOUNT )=  V0RT1 +(Z-A)*( V0RT2-V0RT1 )/(E-A) 

XR(IZ, JCOUNT )=  X 

YR(IZ, JCOUNT )=  Y 

JCOUNT=  JCOUNT+1 

IF  (IP.EQ.O)  10,30 

10  IF  (IFIB.EQ.1 )  GO  TO  20 

CALL  PL0T(X,Y,3) 

20  IP=  1 

GO  TO  100 

30  IF  (IFIB.EQ.1 )  GO  TO  40 

CALL  ?LCT(X,Y,2) 

40  GO  TO  100 

50  XR (IZ , JCOUNT )=  XI 

YR(IZ, JCOUNT )=  Y1 
sxxr(iz, JC0UNT)=  sxxi 
SYYR (IZ, JCOUNT )=  SYY 1 
SXYR(IZ, JCOUNT )=  SXY1 
UR (IZ, JCOUNT )=  U1 
VR(IZ, JCOUNT )=  VI 
V0RTR(IZ, JCOUNT )=  V0RT1 
JC0UNT=  JCOUNT+1 
XR(IZ, JCOUNT)=  X2 
YR(IZ, JCOUNT )=  Y2 
SXXR(IZ, JCOUNT )=  SXX2 
SYYR (IZ, JCOUNT )=  SYY 2 
SXYR(IZ , JCOUNT )=  SXY2 
UR(IZ, JCOUNT)=  U2 
VR(IZ,COUNT)=  V2 
V0RTR(IZ, JCOUNT )=  V0RT2 
JC0UNT=  JCOUNT+1 
IF  (IP.EQ.O)  60,80 

60  IF  (IFIB.EQ.1 )  GO  TO  70 

CALL  PL0T(X1 ,Y1 ,3) 

70  IP-  1 

GO  TO  90 

IF  (IFIB.EQ.1)  GO  TO  90 
CALL  PL0T(X1 ,Y1 ,2) 
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90 

IF  (IFIB.EQ. 1 )  GO  TO 

CALL  PL0T(X2,Y2,2) 

100 

CONTINUE 

RETURN 

END 

Q*********************************************************************** 

c 

c 

c 

c 

SUBROUTINE  ARRAN G ( IC OUN T , NL IN ) 

C 

C 

C 

c 

c********** 


Q********** 

C  ARRANGES  STREAMPOINTS  IN  ORDER  TO  FORM  A  UNIQUE  STREAMLINE 

c********** 


C********** 


c 

c 

c 

c 

c 


DIMENSION  IC0UNT(20) 
COMMON  /A5/SXXR(20,1 ) 
COMMON  /A6/SYYR(20, 1 ) 
COM  ON  /A7/SXYR(20,1  ) 
COMMON  /A9 /UR (20,1) 
COMMON  /A10/VR(20, 1  ) 
COMMON  / A1 1 /YORTR (20,1) 
COMMON  /A12/XR(20, 1 ) 
COMMON  /A1 3/YR(20 , 1 ) 
COMMON  /A14/SXXAR(20, 1 ) 
COMF1  ON  /A1  5/SYYAR(20, 1 ) 
COMMON  /A1 6/SXYAR ( 20 , 1 ) 
COMMON  /A1 7/UAR(20 , 1 ) 
COMMON  /A18/VAR(20, 1  ) 
COMMON  /A 1 9 /VORT AR (20,1) 
COMMON  /A20/XAR(20 , 1 ) 
COMMON  / A2 1 /YAR(20 , 1  ) 
COMMON  /ICOPY/ICOPY , IF IB 


FIND  THE  INITIAL  POINT  IN  THE  STREAMLINE 


DO  90  1=  1 ,NLIN 

DO  10  J=  1 ,IC0UNT(I)-1 

IF  (XR(I, J).NE.O)  GO  TO  10 
XAR(I , 1 )=  XR(I , J) 

YAR(I , 1 )=  YR ( I , J ) 

SXXAR(l , 1 )=  SXXR(I,J) 
SYYAR(I , 1 )=  SYYR(I,J) 
SXYAR(I , 1 )=  SXYR(I, J) 


o  o  o 
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UAR (1,1 )=  UR ( I , J ) 

VAR (I , 1 )=  VR( I , J) 

V0RTAR(I,1)=  V0RTR(I , J) 

GO  TC  20 

1 C  CONTINUE 

C 

C 

C  DETERMINE  THE  NEIGHBOR  OF  THE  INITIAL  POINT 

C 

C 

20  T=  J 

IF  ( T/2 . . EQ . INT ( T/2 . ) )  GO  TO  30 
XAR ( 1 , 2 )  =  XR(I, J+1 ) 

YAR ( I , 2 ) =  YR(I, J+1 ) 

SXXAR(I,2)=  SXXR(I, J+1 ) 

SYYAR(I,2)=  SYYR(I, J+1 ) 

SXYAR(I,2)=  SXYR (I , J+1  ) 

UAR(I,2)=  UR (I , J+1 ) 

VAR (I ,2)=  VR(I, J+1 ) 

V0RTAR(I,2)=  VORTR(I , J+1 ) 

XR ( I ,  J )  =  0. 

YR(I , j)=  0. 

XR(I, J+1 )=  0. 

YR(I, J+1 )=  0. 

C-0  TO  40 

30  XAR ( 1 , 2 ) =  XR(I , J-1  ) 

YAR(I,2)=  YR(I, J-1  ) 

SXXAR (1 , 2 ) =  SXXR(I , J-1  ) 

SYYAR(I,2)=  SYYR (I , J-1 ) 

SXYAR(I,2)=  SXYR (I, J-1  ) 

UAR ( I , 2 ) =  UR(I, J-1 ) 

VAR( 1 , 2 )=  VR(I, J-1  ) 

V0RTAR(I,2)=  VORTR(I , J-1  ) 

XR(I , J)=  0. 

YR  (I ,  J)=  0. 

XR(I, J-1 )=  0. 

YR(I, J-1 )=  0. 

40  CONTINUE 

C 

C 


DETERMINE  THE  REMAINING  STREAMLINE  POINTS  IN  ORDER 


C 

c 

TOL=  .00001 

DO  80  L=  3 > ( IC0UNT(I )-1 )/2+1 
DO  50  M=  1 , IC0UNT(I)-1 

IF  (XR(I,M)+T0L.GE.XAR(I,L-1 ) . AND.XR(l ,M)-T0L. 

LE .XAR (I , L-1 ) . AND .YR(I,M)+TOL.GE. YAR ( I , L- 1 ) . 
AND.YR(I,M)-T0L.LE.YAR(I,L-1  ))  C-0  TO  60 


1 

1 


50 
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CONTINUE 
GO  TO  80 
C 
C 

C  DETERMINE  THE  NEIGHBORS 

C 

C 

60  s=  M 

IF  ( S/2 . . EQ . INT (S/2 . ) )  GO  TO  70 
XAR(I,L)=  XR(I,M+1 ) 

YAR(I,L)=  YR(I ,M+1 ) 

SXXAR(I , L)=  SXXR(I,M+1 ) 

SYYAR(I , L)=  SYYR(l,M+l) 

SXYAR(I , L)=  SXYR (I ,M+1 ) 

UAR(I,L)=  UR(I ,M+1 ) 

VAR(I ,L)=  VR(I ,M+1 ) 

VORTAR(I,L)=  VORTR(l ,M+1 ) 

XR(I ,M)=  0. 

YR(l , J)=  0. 

XR(I ,M+1 )=  0. 

YR(I,M+1 )=  0. 

GO  TO  80 

70  XAR(I,L)=  XR(I,M-1) 

YAR(l , L)=  YR(I,M-1 ) 

SXXAR(I , L)=  SXXR(I.M-l) 

SYYAR(I ,L)=  SYYR(I ,M-1 ) 

SXYAR(I , L)=  SXYR (I ,M-1 ) 

UAR (I , 1)=  UR(I,M-1 ) 

VAR(I,L)=  VR(I,M-1 ) 

VORTAR(l , L)=  V0RTR(I,M-1 ) 

XR(I,M)=  0. 

YR(I,M)=  0. 

XR(I,M-1 )=  0. 

YR(I ,M-1 )=  0. 

C 

C 

80  CONTINUE 

90  CONTINUE 

IF  (IC0PY.EQ.1  )  GO  TO  110 

OPEN (UNI T=3 1 , DEVI CE= ' DSK ’ ,FILE= ' STREAM . DAT ’ ) 

DO  100  J=  1 , NLIN 

VRITE(J1 ,1000) 

DO  100  1=  1 ,(IC0UNT(J)-1 )/2+1 

WRITE(31 ,1010)  XAR ( J , I ) , YAR ( J , I ) , SXXAR ( J , I ) , SYYAR ( J , 
1  I ) ,SXYAR( J ,1) ,UAR( J, I) , VAR ( J, I) , 

1  VORTAR( J, I) 


100 

CONTINUE 

CONTINUE 

CL0SE(UNIT=31 ) 

1 10 

RETURN 

1000 

FORMAT (//////////////////) 
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1010  FORMAT ( 1 X , 8E 1 2 . 4 ) 

ENT 

C 

C 

c 

c 

SUBROUTINE  INTEG ( ICOUNT ,  NUN ) 

C 

c 

c 

c 

0********** 

0********** 

C  INTEGRATES  THE  ORIENTATION  ALONG  A  STREAMLINE 
0********** 

0********** 

DIMENSION  IC0UNT(20) 

DIMENSION  VAL(20) 

DIMENSION  VAT (20) 

DIMENSION  Y(1 ),C(24),W(1  ,9) 

COMMON  /A14/SXXAR(20,1 ) 

COMMON  /A 1 5/SYYAR ( 20 , 1 ) 

COMMON  /A1 6/SXYAR(20, 1 ) 

COMMON  /A1 7/UAR (20 , 1 ) 

COMMON  / A1 S/VAR (20 , 1 ) 

COMMON  /A1 9/VORTAR (20,1) 

COMMON  /A20 /XAR (20,1) 

COMMON  /A21 /YAR ( 20 , 1 ) 

COMMON  /A22/I , J, B 
COMMON  / A26/PHI (10,20,60) 

COMMON  /ICOPY /ICOPY, IF IB 
COMMON  /ITER/ITER, OCOUNT 
COMMON  /IHERM/IHERM 
EXTERNAL  FCN1 

DATA  (VAL(I ) , 1=1 , 20 ) /— 1 .50,-1 -50,-1 .50,-1 .50,-1 .50,-1 .50,-1.50, 
1  -i .50,-1 -50,-1 .50,-1 .5,-1 -5,-1 -5,-1 -5,-1 -5,-1 -5,-1 -5,-1 -5, 

1  -1.5, -1-5/ 

DO  10  1=  1,20 

VAT(I)=  VAL(I )  +  (OCOUNT-1 )/ITER*5- HI  592654 
10  CONTINUE 

IP  ( ICOPY. EQ.1)  GO  TO  20 
20  CONTINUE 

B=  1  .0 

DO  60  1=  1.NLIN 
IND=  1 
TOL=  .0001 
NW=  1 
N=  1 

Y(1 )=  V AT ( I ) 

XEND=  0.0 


X=  0.0 

DO  50  J=  1 ,(IC0UNT(I)-1 )/2+1 
IF  ( J.EQ. 1 )  GO  TO  50 

TOP=  (( (XAR(I, j)-XAR(I, J-1 ) )**2 . +  (YAR ( I , J)-YAR(I , J- 
1  1 ))**2.)**«5) 

BOTTOM=  ( (UAR(I , J-1 )**2+VAR(I , J-1 )**2)**.5) 

TINC=  TOP/BOTTOM 
XEND=  XEND+TINC 

CALL  DVERK(N,FCN1 , X , X , XEND , TOL , INB , C , NW , W , IER ) 

50  IF  (IHERM.EQ. 1 )  GO  TO  40 

CALL  FPLOT(I,J,Y) 

40  CONTINUE 

PHI(0C0UNT, I, J)=  Y(1  ) 

IF  (ICOPY.EQ. 1 )  GO  TO  50 

WRITE(59 ,1010)  Y ( 1 ) , XEND , XAR ( I ,  J ) , YAR ( I ,  J ) , UAR ( I , J ) , 

1  VAR ( I, J),TINC, TOP, BOTTOM 

50  CONTINUE 

IF  (ICOPY.EQ. 1)  GO  TO  60 
WRITE(59 , 1 020) 

60  CONTINUE 

IF  (ICOPY.EQ. 1 )  GO  TO  70 
70  RETURN 

1000  FORMAT (5X, II ) 

1010  FORMAT ( 1 X , 9E 1 0 . 5 ) 

1020  FORMAT(// /////////////) 

END 

Q**************************************************************** ******* 

c 

c 

c 

c 

SUBEOUTINE  FCN1 (N,X, Y, YPRIME) 

C 

C 

c 

c 

0**** ****** 

0********** 

C  EXTERNAL  SUBROUTINE  FOR  INTEG 
0********** 

0********** 

DIMENSION  Y ( 1 ) , YPRIME ( 1  ) 

COMMON  /A1 4/SXXAR (20 , 1 ) 

COMMON  /A15/SYYAR(20,1 ) 

COMMON  /A16/SXYAR(20,1 ) 

COMMON  /A17/UAR(20,1 ) 

COMMON  /A1 8/VAR (20,1 ) 

COMMON  /A19/V0RTAR(20,1 ) 

COMMON  /A20 /XAR (20,1 ) 

COMMON  /A21 /YAR (20,1 ) 

COMMON  /A22/I , J,B 
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10 

1000 


COMMON  /ICOPY/ICOPY , IF  IB 

YPRIME ( 1 )=  VORTAR( I , J)+B*(SXYAE ( I , J 

1  ( SYYAR (I , J) /2 . DC-SXXAE ( I 

IF  (ICOPY.EQ. 1 )  GO  TO  10 

RETURN 

FORMAT (20X.4E1 0.3) 

END 


/2. EC*CGS (2*Y ( 1 ) ) -0 . 5* 
J)/2 • LC)*SIN(2*Y ( 1 ))) 


c 

c 

c 

c 

SUBROUTINE  BGUNP(ICOPY , SC ) 
C 
C 
C 


C 

0********** 

0********** 

C  SUBROUTINE  TO  PLOT  THE  BOUNDARY 

C********** 

Q********** 

COMMON  / A1 /NNCDE , NVAR , NELEM , NBD , NSO , NVEL 
COMMON  /A2/X ( 1 ) 

COMMON  /A3/Y (1  ) 

COMMON  /A23/EX , BY , MINX , MINY , XMIN , YMIN , SCALE , MENU 
C  $  ***INTRODUCE  MODE  TWICE*** 

Q********4('******4i'*************************************  ******** 

DIMENSION  XP(254) ,YP(254) 

DO  10  1=  1 , NNODE 

XP(I)=  BX*X(I) 

YP(I)=  BY*Y ( I ) 

10  CONTINUE 

DO  20  1=  1 , NNODE 

XP(I)=  MINX+(XP(I )-XMIN )*SCALE 
YP(I)=  MINY+(YP(I)-YMIN)*SCALE+1 . 

20  CONTINUE 

IF  (ICOPY.EQ. 1)  30,40 
30  CALL  PLTSRTC JONES' ,30) 

GO  TO  50 

40  CALL  PLTSRTC  JONES' ,20) 

50  CALL  FACTOR (SC) 

60  CALL  PL0T(0. ,0. ,3) 

CALL  PL0T(8-5,0. ,2) 

CALL  PL0T(8.5,11 .5,2) 

CALL  PL0T(0. ,1 1 .5,2) 

CALL  PL0T(0. ,0. ,2) 

CALL  PL0T(0. ,1 . ,3) 

CALL  PL0T(8. 5 , 1 . , 2) 

CALL  SYMB0L(0. 5, 0.375, 0.225, MENU, 0. ,30) 
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C 

C  PLOT  OF  THE  BOUNDARY 

C 

C 

XI =  XP(1 ) 

II =  YP(  1 ) 

CALL  PLOT(X1 ,Y1 ,3) 

NBD2=  NBD/2 
DO  70  1=  1 , NBD2 
J=  1  +  1 

IF  (I.EQ.NBD2)  J=  1 
XXX=  XP(J) 

YYY=  YP( J) 

70  CALL  PLOT (XXX , YYY , 2 ) 

CONTINUE 

RETURN 

END 

Q*********************************************************************** 

c 

c 

c 

c 

SUBROUTINE  FPLOT(l,J,Y) 

C 

c 

c 

c 

0********** 

0********** 

C  PLOTS  INDIVIDUAL  FIBER  ORIENTATIONS 

0********** 

0********** 

DIMENSION  XART ( 20 , 1 00 ) , YART (20,100) 

DIMENSION  Y(1) 

COMMON  /A23/BX , BY , MINX ,MINY , XMIN , YMIN , SCALE 
COMMON  /A20/XAR(20 , 1 ) 

COMMON  /A21 /YAR(20 , 1 ) 

COMMON  /A30/HEIGHT 
COMMON  /IHERM/IHERM 
IF  (IHERM.EQ. 1 )  GO  TO  10 
PI=  3.U1592654 
FLENG=  0.025*HEIGHT 
XART(I, J)=  XAR(I.J) 

YART(I , J)=  YAR(I,J) 

C 

C 

C  RESCALE  TO  PLOTTING  COORDINATES 

C 

C 
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XART ( I ,  J )=  MIEX+(XART(I ,  j)*£X-X!4IK)*SCALE 
YART (I , J)=  K INY + ( YART ( I , J ) *EY- YMIU )*SCALE h 


II 

Ph 

X 

XART 

d,j 

)  -FLEUG 

*COS 

(PI 

/2 • -Y ( 1 

J  J 

II 

YART 

(i,j 

)-FLELG 

*SIN 

(PI 

/2 .  —  Y  ( 1 

)  ) 

CALI 

,  FLO 

T(XF 

,YP,3) 

XL  = 

XART 

(I,J 

) +FLENG 

*COS 

(PI 

/2. -Y( 1 

)) 

YL= 

YART 

(I,J 

)+fleng 

*SIN 

(PI 

/2 .  —  Y  ( 1 

)) 

CALL  PLOT(XL , YL , 2 ) 

10  RETURN 

END 

q******************************************** *************************** 


C 

C 

SUBROUTINE  VALKER ( ICCUNT , NLIN ) 
C 
C 
C 


0********** 

0********** 

C  CALCULATES  VALUES  OF  HERMANS  ORIENTATION  PARAMETER 
0********** 

0********** 

DIMENSION  IC0UNT(2C) 

COMMON  /A26 /PHI (10, 20, 60) 

COMMON  /ICOPY/ICOPY , IF  IB 
COMMON  /A27/PHIL0C (20,60) 

COMMON  /ITER/ ITER, OCOUNT 
COMMON  /A28/HERM(20,60) 

COMMON  /A20/XAR(20, 1  ) 

COMMON  /A21/YAR(20,1 ) 

IF  (ICOPY.EQ. 1 )  GO  TO  10 

OPEN (UNIT=47 , DEVICE= ' DSK 1 , FILE =' HERM.DAT’ ) 

10  DO  50  N=  1 , NLIN 

BO  40  J=  1 ,(IC0UNT(N)-1 )/2+1 
PHIM=  0.0 
DO  20  1=  I, ITER 

PHIM=  PHIM+PHI ( I , N , J) 

20  CONTINUE 

PHIM=  PHIM/ITER 
PHSQ=  0.0 
DO  50  1=  1 , ITER 

PHSQ=  PHSQ+(PHI(I,N,J)-PHIM)**2. 

50  CONTINUE 

PHIRMS=  (PHSQ/ITER)** . 5 
PHILOC(N,J)=  5. 141 592654/2 -PHIM 
HERM(N , J)=  2*((C0S(PHIRMS))**2)-1 . 

IF  (HERM(N,J).LT.O)  PHILOC(N,J)=  5. HI  592654/2. 

-PHILOC(N, J) 


1 


197 


IF  (HERM(N, J) .LT.C)  KERK(N,J)=  -KEEM ( N ,  J ) 

IF  (IC0PY.EQ.1)  GO  TO  40 

WRITE (47, 1000)  XAR(N, J) , YAR(N, J) ,HERM(N , J) , PHILOC (N . 

1  J ) , FHIRMS , PHIM , PHSQ 

40  CONTINUE 

50  CONTINUE 

IF  (ICOPY.EQ. 1 )  GO  TO  60 
CL0SE(UNIT=47 ) 

60  RETURN 

1000  FORMAT (7E 14. 6) 

END 

Q*********************************************************************** 

c 

c 

c 

c 

SUBROUTINE  HERMP ( NLIN , IC OUNT , HIN , HINCR , NHLIN ) 

C 

c 

c 

c 

Q ********** 

c********** 

C  SCANS  INDIVIDUAL  TRIANGLES  FOR  CONTOURS  OF  HERMANS 
C  ORIENTATION  PARAMETER 

C********** 

^ ********** 

DIMENSION  I COUNT (20) 

DIMENSION  XXXX(5) ,YYYY(5) ,HHH(3) 

COM  ON  /A20/XAR(20, 1  ) 

COMMON  /A21/YAR(20,1 ) 

COMMON  /A28/HERM(20 , 60) 

COMMON  /ICOPY/ICOPY, IF IB 
DO  70  HIZ=  1, NHLIN 

HLIN=  HIN+HINCR*(HIZ-1  ) 

IF  (ICOPY.EQ. 1)  GO  TO  10 
WRITE (5, 1000)  HLIN 
10  CONTINUE 

DO  60  1=  1 , NLIN-1 

MIN=  (ICOUNT(I)-I )/2+1 

IF  (IC0UNT(I+1 ) .LE.ICOUNT(I))  MIN=  (IC0UNT(I+1 )-1 )/2+ 

1  1 
DO  20  J=  1 ,MIN-1 

XXXX(1 )=  XAR(I, J) 

YYYY(1 )=  YAR(I, J) 

HHH(1)=  HERM(I, J) 

XXXX(2)=  XAR(I+1 , J) 

YYYY(2)=  YAR(I+1 , J) 

HHH(2)=  HERM(I+1 ,J) 

XXXX(3)=  XAR(I , J+1 ) 

YYYY(3)=  YAR(I, J+1 ) 


198 


HHH( 3  )  =  HERK(I , J+1 


Li,  Xi 


.  Ai.h  ( HLlii ,  Xaaa 


.  YYYY . 


XXXX(1 )=  XAR (1+1 , J) 

YYYY(1 )=  YAR (1+1 , J ) 

HKK( 1 )=  HERM (1+1 , J ) 

XXXX ( 2 ) =  XAR (1+1 , J+1 ) 

YYYY(2)=  YAR ( 1+1 , J+1 ) 

HHH ( 2 ) =  HERM(I  +  1 , J+1  ) 

XXXX ( 3 )  =  XAR (I, J+1 ) 

YYYY(3)=  YAR(I , J+1  ) 

KHK ( 3 )  =  HERM(I , J+1 ) 

CALL  TRIANH ( KLIK , XXXX , YYYY , KKH ) 

20  CONTINUE 

IF  (IC0UNT(I+1 ) .EQ.ICOUNT(I))  GO  TO  60 
IF  ( IC OUi'IT ( I  + 1  ) .  LT .  ICOUNT(I ) )  GO  TO  40 
DO  30  H=  J, (IC0UNT(I+1 )-1  )/2 
XXXX(1)=  XAR ( I , J ) 

YYYY(1 )=  YAR ( I , J ) 

HHK( 1 )=  KERM ( I , J ) 

XXXX(2)=  XAR (1+1 ,M+1 ) 

YYYY(2)=  YAR (1  +  1 ,M+1  ) 

HHK(2 )=  HERM ( 1  +  1 ,M+1  ) 

XXXX(3)=  XAR  (1  +  1 
YYYY(3)=  YAR(I+1 ,M) 

HHH ( 3 )  =  HERM(I  +  1,M) 

CALL  TRIANH(HLIN , XXXX , YYYY , HHH ) 


30 

CONTINUE 

C-0  TO  60 

40 

DO  50  M=  J, (ICOUNT(I)-I )/2 

XXXX(1)=  XAR ( I , M ) 

YYYY(1 )=  YAR (I ,M) 

HHH(1  )=  HERM( I ,M) 

XXXX(2 )=  XAR ( I , M+ 1 ) 

YYYY(2)=  YAR (I ,M+1 ) 

HHH ( 2 ) =  HERM(I ,M+1 ) 

XXXX(3)=  XAR (1+1 , J) 

YYYY ( 3 ) =  YAR (1+1 ,J) 

HHH(3)=  HERM(I+1 , J) 

CALL  TRIANK(HLIN, XXXX, YYYY, HHH) 


50 

CONTINUE 

60 

CONTINUE 

70 

CONTINUE 

RETURN 

1000 

FORMATC  CONTOUR=' ,E9-3) 

END 

c 

c 

c 

c 
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SUBROUTINE  TRIANH (HLIN , XXXX , XYYY , HHH ) 

C 

C 

C 

C 

c********** 

c********** 

C  PLOTS  CONTOURS  OF  HERMANS  ORIENTATION  PARAMETER 

Q ********** 

Q ********** 

DIMENSION  XXXX(3),YYYY(3),HHH(3) 

COMMON  /A23/BX , BY , MINX , MINY , XMIN , YMIN , SCALE , MENU 
IP*  0 

DO  10  1=  1 ,3 
C 
C 

C  RESCALE  TO  PLOTTING  COORDINATES 

C 

C 

XXXX(I)=  MINX+(XXXX(I )*BX-XMIN)*SCALE 
YYYY(I )=  MINY+(YYYY(I)*BY-YMIN )*SCALE+1 . 

10  CONTINUE 

DO  80  M=  1  ,3 
K1=  M 

K2=  M+1 -3*(M/3) 

A=  HHH(K1 ) 

B=  HHH(K2) 

C=  (HLIN-A)*(HLIN-B) 

IF  (C.GT.O)  GO  TO  80 
XI =  XXXX(K1 ) 

X2=  XXXX(K2) 

Y1 =  YYYY(K1 ) 

Y2=  YYYY(K2) 

IF  (A.EQ.B)  GO  TO  40 

X=  XI +(HLIN-A)*(X2-X1 )/(B-A) 

Y=  Y1 +(HLIN-A)*(Y2-Y1 )/(B-A) 


20 

IF  (IP.EQ.O)  20,30 
CALL  PL0T(X,Y,3) 

30 

IP*  1 

GO  TO  80 

CALL  PL0T(X,Y,2) 

40 

GO  TO  80 

IF  (IP.EQ.O)  50,60 

50 

CALL  PLOT (XI  ,Y1 ,3) 

60 

IP=  1 

GO  TO  70 

CALL  PLOT (XI ,Y1 ,2) 

70 

CALL  PL0T(X2,Y2,2) 

80 

CONTINUE 

RETURN 

END 

